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Preface 


Methods for numerical simulations in engineering sciences keep updating every 
day, and naturally the state-of-the-art advancement may not be incorporated into the 
well-known textbooks. This short book introduces the Conservation Element and 
Solution Element (CESE) method, which has a relatively short history. The CESE 
method provides time-accurate simulations of physical systems and overcomes the 
difficulties in capturing discontinuities, by combining many innovative considera- 
tions in a nontraditional manner. This book elaborates on how the CESE method 
works and presents a comprehensive survey of the recent CESE developments and 
applications. 
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Chapter 1 A) 
Introduction E 


With the rapid development of electronic computers, numerical computation has 
become an important paradigm of scientific discovery as well as a powerful tool for 
engineering research. Solving complex problems in a computational fashion is more 
than applying theories, equations, and formulas. Computational methods, also called 
algorithms or schemes, have strong influences on the outcomes of computations. 

The space-time conservation element and solution element (CESE) method is 
aimed to numerically solve equations of conservation laws in various physical 
systems. The major concern herein will be the computational fluid dynamics (CFD), 
which is a hot area in scientific and engineering researches. Nevertheless, it may help 
at the outset to recognize that the conservation laws in CFD problems have a lot in 
common with conservation laws in other fields, such as acoustics, solid dynamics, 
electromagnetics, and magnetohydrodynamics. The physics may be different, but 
the mathematics are similar. For instance, physics involving dynamical evolution 
of waves and discontinuities are usually modelled by time-dependent nonlinear 
hyperbolic partial differential equations (PDEs). Some CFD problems happen to 
be representative of such physical problems in a wider context. 


1.1 Background 


Most of the key achievements in conventional CFD have been incorporated into the 
computational procedure of the finite volume method (FVM). Indeed, the FVM is 
well established and widely used, and therefore it is sometimes regarded as a mature 
technique. There are two important steps in the FVM, namely the reconstruction and 
the evolution [1, 2]. In the reconstruction step, one needs to locally approximate the 
flow field with simple functions such as polynomials, and then interpolate the values 
of flow quantities at the cell interfaces by utilizing the cell-average values. In the 
evolution step, the numerical flux at each cell interface is determined by some kind 
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2 1 Introduction 


of flux technique, and then the time integration of the conservation equations can be 
performed. 

A major category of CFD problems face the challenges brought by (1) the coexis- 
tence of smooth regions and discontinuities in the flow (e.g., shock waves and material 
interfaces) and (2) the wide spectrum of wave numbers/frequencies in the flow (e.g., 
acoustic waves and turbulence). Due to these physical phenomena, the CFD solver is 
required to well handle the numerical dissipation, which leads to an endless debate in 
the CFD community. First, for the robustness in capturing shock waves, there must 
be a certain amount of numerical dissipation to suppress the spurious oscillations. 
However, a low numerical dissipation is very favourable for the accuracy and reso- 
lution of the viscous layers, contact surfaces, and small-scale structures in the flow, 
which would be smeared out by excessive dissipation. 

To achieve the low but controllable numerical dissipation, the research works 
on FVM are mostly dedicated to two approaches: (1) developing high-order recon- 
struction techniques in conjunction with nonlinear limiters and (2) improving the 
flux solvers (e.g., approximate Riemann solvers). In both ways, numerous important 
advances have been obtained. 

In spite of these successes in the FVM development, several shortcomings may 
arise. (1) A conventional FVM implementation highly relies on the selections of 
special techniques, such as the limiter and the Riemann solver. Actually, each special 
technique has its own pros and cons. There are many alternatives for each option, 
but none of them can be proven to be optimal and general. (2) The coupling of space 
and time is usually overlooked and may not be guaranteed in the numerical solution. 
(3) The multi-dimensionality can be questionable because most flux solvers are only 
based on one-dimensional physics. (4) The popular high-order FVM schemes are 
usually not compact, i.e., large stencils are used. Under this background, the CESE 
method was proposed to overcome the difficulties in the conventional FVM to some 
extent, and it became a new member in the family of FVM. 


1.2 History of the CESE Method 


The space-time conservation element and solution element method was initiated 
by Sin-Chung Chang of NASA Glenn Research Centre and his collaborators for 
aerodynamic computation in early 1990s [3, 4]. There were several motivations to 
propose the CESE method. First, they wanted to construct multidimensional and 
space-time unified CFD schemes, without using dimensional splitting or separated 
treatments of space discretization and time discretization. Second, they believed a 
CESE solver should be built from a non-dissipative core scheme so that numerical 
dissipation can be controlled effectively, dynamically, and even actively. Third, they 
attempted to avoid using the Riemann solver in the CESE method. 

Around 1994, the CESE project was approved and supported by NASA Glenn 
Research Centre. The work of the research group showed that this time-accurate 
scheme possesses low numerical dissipation, which is valuable in CFD simulation. 
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For the first time, the transonic resonance in a convergent-divergent nozzle with 
frequency-staging effect was successfully simulated by the CESE scheme [5, 6]. In 
this case, the blind prediction by the CESE scheme matched the experimental data. 
Therefore, the CESE method is proved to be capable of simultaneously simulating 
multidimensional unsteady shock waves and acoustic waves with high accuracy. 

Owing to its accuracy and robustness, the CESE method is employed in the 
CFD module of the simulation software LS-DYNA [7]. The Jacobs Technology Inc. 
(designer of the NASA’s hypersonic flow test facilities) developed its own in-house 
CESE code called JUSTUS (Jacobs Unified Space-Time Unstructured Solver) for 
hypersonic flow simulations. Another in-house CESE code, called “ez4d” software, 
was developed by the NASA Langley Research Centre. This is a time-accurate 3D 
Navier-Stokes flow solver on unstructured meshes [8, 8]. 


1.3 Main Features of the CESE Method 


The CESE method is a special finite-volume-type numerical method, with the aim 
of solving the governing equations of fluid dynamics and other conservation laws in 
various physical systems. The CESE method possesses many features such as the 
unified treatment for space and time, the fully discrete one-step explicit scheme, and 
the highly compact stencil. Without enlarging the stencil or adding stages of time 
integration, the CESE method can achieve arbitrary high-order accuracy for both 
space and time. 

The essential ingredients in the CESE method include: (1) the spatial deriva- 
tives of physical variables are stored as independent unknowns, in addition to the 
physical variables themselves. In every time step, these derivatives are updated by a 
specially designed procedure. (2) a staggered mesh and a staggered time-marching 
strategy are employed. (3) the interior structure within each solution element is 
built with the Taylor expansion. (4) the time-marching approach is based on the 
Cauchy—Kowalewski procedure. 

Through a combination of the above techniques, the CESE method has low dissi- 
pation and high compactness. When applied to the simulations of complex physical 
processes, the CESE method can catch shock waves, contact discontinuities, fine 
structures, and small disturbances, with high resolution and strong robustness. There- 
fore, the CESE method demonstrates good performances in the numerical simula- 
tions of wave-propagation problems (e.g., shock waves, acoustic waves, detonation 
waves, stress waves in solid, and the electromagnetic waves), interfacial instabil- 
ities, as well as the interactions of gaseous and liquid phases. In many research 
areas including high-speed aerodynamics, shock dynamics, detonation, aeroacous- 
tics, and solid dynamics, the CESE method proves to be suitable and shows a good 
development prospect. 


4 1 Introduction 


1.4 Outline of the Book 


The remainder of this book is organized as follows. First of all, the non-dissipative 
core scheme of the CESE method is introduced in Chap. 2, with detailed descrip- 
tions of the basic concepts. Then, the practical shock-capturing CESE schemes with 
numerical dissipation including the classical a—a scheme, the Courant number insen- 
sitive scheme, and the recently proposed upwind CESE schemes are presented in 
Chap. 3. Furthermore, Chaps. 4 and 5 extend the CESE method to multidimensional 
and high-order versions. In Chap. 6, numerical properties of various CESE schemes 
are analysed, along with comparisons to other numerical schemes. Chapters 7—9 
provide an overview of the applications of the CESE method, most of which are 
quite relevant to the aerospace engineering. Finally, a summary and an outlook of 
the CESE method are given in Chap. 10. 
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Chapter 2 R) 
Non-dissipative Core Scheme of CESE checi; 
Method 


This chapter is devoted to demonstrating the basic ideas in the CESE method. These 
ideas include the adoption of a space-time integral form of governing equations as 
the starting point of scheme construction, as well as the introduction of conservation 
element (CE) and solution element (SE) in the discretization of space-time domain. 
Then, the non-dissipative core scheme of the CESE method will be presented in 
detail. 


2.1 Space-Time Integral Form of Governing Equations 


Governing equations for a specific physical problem can be expressed in both differ- 
ential and integral forms. Usually, a differential equation or a set of differential 
equations can clearly describe the evolution of a physical system. The differential 
forms of governing equations are prevalent in textbooks and theoretical research 
works, due to their compactness in writing and the maturity in the mathematical 
analysis on differential equations. 

In the area of numerical simulation, the finite difference method is directly applied 
to the differential equations, while the finite volume method and finite element 
method require the governing equations to be cast into integral forms before numer- 
ical treatments. In fact, the performance of a numerical method will be influenced 
by the form of equations that is used in the numerical method, even though different 
forms of governing equations can be mathematically equivalent. As will be shown 
in this book, a major feature of the CESE method is the adoption of a space-time 
integral form of the governing equations, in which time and space are treated on the 
same footing. 

For illustrative purposes, we consider the compressible Euler equations for 2D 
planar flows, which can be written in a differential form as 
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aU əF ðG 


—+—+4—=0, 2.1 
ot ox dy on 


where U denotes the vector of conserved variables such that 


U= ; (2.2) 


where F and G are the inviscid fluxes of the form 


pu pv 
2 
puv 
pa | TP | ge ; (2.3) 
puv pv +p 
(E + p)u (E+ p)v 


In Eqs. (2.2) and (2.3), p is the density of the fluid, u and v are the x-component 
and y-component of the flow velocity, respectively, p is the static pressure, and E is 
the total energy per unit volume of the fluid. Note that once the equation of state 
(EOS) is provided, Eq. (2.1) becomes a closed set of equations. Here, flux vectors 
F and G are regarded as functions of conserved vector U. A set of equations in the 
form of Eq. (2.1) is called a system of conservation laws. 

By using the gradient operator in 2D physical space 


a ə 
V= (è x). (2.4) 


Equation (2.1) can also be written in a divergence form 


aU 
V-H= 2. 
z + 0, (2.5) 


where H is a matrix composed of spatial flux vectors F and G: 
H = (F, G) (2.6) 


Therefore, the divergence form of Euler equation states that at time t and point 
(x, y), the temporal rate of change of conserved quantities plus the divergence of the 
spatial flux of these quantities must be zero. 

Consider a control volume denoted by V in 2D physical space, as shown in Fig. 2.1, 
with the surface S(V) and the unit outward normal vector n on the surface. Integrating 
Eq. (2.5) over the control volume V and applying the Gauss’s theorem, one obtains 
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Fig. 2.1 Schematic of a 
control volume in 
two-dimensional physical 
space 


O 


ə 
+([ wav) + H-ndS=0, (2.7) 
ðt \Jy SV) 


which can be interpreted as follows: the temporal rate of change of conserved quan- 
tities within the volume V must equal to the net inward flow rate of these quantities 
through the surface S(V). This integral form of Euler Eq. (2.7) is the starting point 
of the well-known finite volume method (FVM) in computational fluid dynamics. 

With a unified treatment of time and physical space, we can define x, y, and t as 
coordinates in a three-dimensional Euclidean space, called the space-time domain. 
On this basis, the definition of gradient operator can be extended as 


ð ə ə 
V= (ż =, ) (2.8) 
əx oy ot 


and Eq. (2.1) can be written in a divergence-free form as 
V-h=0, (2.9) 
where the matrix h is composed of both flux vectors and the conserved vector U: 
h = (F, G, U) (2.10) 
By applying Gauss’s theorem to the integration of Eq. (2.9) over an arbitrary 


control volume V in the 3D space-time domain as shown in Fig. 2.2, the Euler 
equation is eventually formulated as 


§ h-nds =0 (2.11) 
S(V) 


which means the balancing of flux is enforced for a space-time control volume. 
Clearly, this equation gives a direct description of the space-time conservation of 
mass, momentum, and energy in fluid flows, which faithfully preserves the original 
physical laws. 
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Fig. 2.2 Schematic of a Í 
control volume in 
three-dimensional 
space-time domain 


X 


Equation (2.11) provides an example of the space-time integral form of governing 
equations for physical problems. It is worth noting that the space-time integral form 
and the traditional integral form like Eq. (2.7) state the same physics and must be 
equivalent to each other in mathematics, but they can lead to different numerical 
schemes. In contrast to the finite volume method, the CESE method numerically 
implements the space-time integral form of governing equations, emphasizing on 
the conservative property in the unity of time and physical space. 


2.2 Definitions of Conservation Elements and Solution 
Elements 


The CESE method starts with discretizing the space-time domain which is relevant 
to the computation. The discretization procedure generates the mesh for the physical 
space and determines the time-step size for the time-marching algorithm, similar to 
most CFD methods. The features of the CESE discretization include the construction 
of control volumes for the space-time integral form of the governing equation, the 
arrangement of the solution points, and the selection of the unknown variables to be 
calculated and stored at each solution point. All these features can be demonstrated 
by introducing two special concepts: the conservation element (CE) and the solution 
element (SE), which coin the algorithm accordingly. 

To explain the CESE method in a simple way, we begin with the application of the 
CESE method to a 1D scalar conservation law of the form, without loss of generality 


ot Ox 


“poi (2.12) 


With a uniform division of the 1D physical space and a constant time step, the 
space-time mesh is constructed as a 2D x-t plane as shown in Fig. 2.3 (solid lines). 
The spatial coordinate of the j-th mesh node is denoted by x;, with the corresponding 
cell centre position at Xj41/2 = (x; + xj41)/2, and the cell size of Ax = xj4, — x;. In this 
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| x 


Fig. 2.3 Mesh and arrangement of solution points for 1D CESE scheme 


time-marching algorithm, each step between źt„,—ı and ¢, consists of two half steps: 
[tn—1, tn—12] and [tp 1/2, tn], with the time-step size of At = tn — ty-1. 

The unknown function u(x, t) is represented by the discrete values of u at a set 
of specific space-time points, called solution points which are shown as mesh nodes 
in Fig. 2.3—black circles for integer time levels {to, t1, ..., fn, ... } and red squares 
(cell centres) for half-integer time levels {t1/2, t3/2, ..., tn+1/2,---}. In other words, a 
staggered mesh is used at each intermediate time level. In the CESE scheme, both 
the unknown variable u and its spatial derivative uy will be calculated and stored at 
each solution point (e.g. point (j, n)): 


ðu 
u; = u(xj, tn), (Ux); = ay tn) (2.13) 


The paths of information flow in a single CESE time step is illustrated in the 
schematic diagram in Fig. 2.4. As seen, a highly compact stencil in space and time 
is used in the CESE scheme. If a half time step it is treated as the basic iteration, the 
halfwidth of the symmetric stencil is then Ax/2 and unknowns at point (j, n) only 
depend on the data stored at (j — 1/2, n — 1/2) and (j + 1/2, n — 1/2). 

Based on the space-time mesh shown in Fig. 2.3, a set of small space-time 
elements, named conservation elements (CE), can be constructed, as demonstrated 
in Fig. 2.5. A CE is assigned to each solution point. For example, (CE); is denoted 
as the CE for point ( j, n), which is the rectangle with four vertices at the points 
G — 1/2, n — 1/2), (j + 1/2, n — 1/2), G + 1/2, n), and (j — 1/2, n). Apparently, 
the CEs cover the whole space-time domain without overlap. It is noteworthy that 
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Fig. 2.4 The information flow in one CESE time step with time-marching variables 


the arrangement of CEs is staggered in two successive half time steps. Accordingly, 
the space-time integral form of Eq. (2.12) is numerically implemented within each 
conservation element and discrete equations for unknowns are established. 

When performing the space-time integration of Eq. (2.12) over a CE, an important 
question arises: how to evaluate u and f along the boundary of the CE. This leads to the 


n+l 


Fig. 2.5 Schematic of conservation elements (CEs) arrangement 
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Fig. 2.6 Schematic of the 
solution element (SE) 
associated with solution 


point (j, n) 


introduction of a solution element (SE) for each solution point, as shown in Fig. 2.6. 
(SEY; composed of two line segments bisecting each other at point (j, n), forming 
a cross with four endpoints (j + 1/2, n), (j, n + 1/2), G — 1/2, n), and (j, n — 1/2). 
For half-integer points, SEs can be defined in the same way. Note that conservation 
element (CE); is bounded by line segments belonging to three solution elements: 


(SE)", (SE); 1/2» and (SE)';, 1/2- In the CESE scheme, the functions u(x, f) and f (x, 
t) are assumed to be linear along each line segment of a SE and are approximated by 
first-order Taylor expansions about the centre of the SE. Specifically, in the solution 


element (SE)’, u and f are constructed as 


u(x,t) =u" + (u(x — xj) + (UCE — th), Œ, t) € (SE)? (2.14) 


Ft) = rUe- + FN = th), (0,1) € (SE) (2.15) 


where u; and f, are the temporal derivatives of u and f, respectively. 


2.3 Non-dissipative Core Scheme: a Scheme 


In this section, we present a non-dissipative CESE scheme, named the a scheme, to 
solve the 1D scalar conservation law, Eq. (2.12). A space-time flux vector is defined 
as 


h= (f,u) (2.16) 


where f and u can be regarded as the components of flux vector h in x-direction and 
t-direction, respectively. According to this definition, Eq. (2.12) can be converted into 
the space-time integral form of Eq. (2.11) by following the procedure in Sect. 2.1. 
Consider a half step marching from time level n— 1/2 to time level n. At the solution 
point (j, n), two independent unknowns uw’; and (ux), need to be calculated simul- 
taneously. Therefore, two algebraic equations need to be formulated by discretizing 
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Fig. 2.7 Definition of (SE j; 
sub-CEs: CE~ and CE*, and J 
the associated SEs 


the integral conservation law. For this purpose, the conservation element (CEY; is 
split into two sub-elements: (CE`); and (CE®);. As shown in Fig. 2.7, (CE )’ and 
(CE*)5 are the rectangle ACDE and the rectangle CBFD, respectively. Each edge of 
these rectangles belongs to one of the three SEs associated with (CE)’. 

Next, the space-time integral conservation law is implemented over each of the 
sub-CEs. Let the control volume V in Eq. (2.11) be (CE )i and (CE+); in turn, one 
can obtain 


$ h -ndS = 0 (2.17) 
S(CE-); 
and 


f h-nds =0, (2.18) 
S(CEt)" 


where h is the space-time flux vector expressed in Eq. (2.16) and n is the unit outward 
normal vector on the boundary of (CE )’ or (CE+);. 

As depicted in Fig. 2.8, the average values of u(x, t) on line segments DE, DF, 
AC, and BC are denoted by U L, Up, Uz, and Up, respectively, while the average 
values of f (x, t) on line segments AE, BF, and CD are presented by Fz, Fr, and Fc, 
respectively. Notably, CD is the interface between two sub-CEs. With this notation, 
integral Eqs. (2.17) and (2.18) can be expressed as 


Ge ai aS (2.19) 
L3 Mi L aos : 

x Ax Ax At 

URS = Un + (Fe FS (2.20) 


which explicitly state the balance of space-time flux in each sub-CE. 
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I 
1 
Fm (CE); 


n 


U; 


Fig. 2.8 Space-time flux through the each boundary of sub-CEs 


To proceed with the scheme construction, the time-marching variables (i.e. u and 
ux) at the solution points will be linked with UL“, Ur’, UL, Ur, Fr, Fr, and Fc 
in Eqs. (2.19) and (2.20), using the concept of the solution element. Recall that u 
and f can be calculated by first-order Taylor expansions about the solution point, by 
assuming them to be linear functions of x and t inside each individual SE. Since line 
segments AC and AE belong to (SE); o U, and Fz can be obtained in terms of the 
known information stored at (j — 1/2, n — 1/2), i.e., point A in Fig. 2.8. Performing 


Taylor expansions in (SEZ yields 


n— Ax PE 

U, = uji + T UV (2.21) 
n—1/2 At n—1/2 

FL = fib + g h (2.22) 


Similarly, because BC and BF belong to (SE) 12 Up and Fp can be obtained 
by the Taylor expansions about solution point (j + 1/2, n — 1/2): 


n= Ax n— 

Ur = wii — qd (2.23) 
n—1/2 At n—1/2 

Fr = Fak + F an (2.24) 


In Egs. (2.21)-(2.24), u and ux at time level n — 1/2 are the known values. In 
addition, since f = f (u) is a prescribed function in the conservation law, 


fa, ETE, (2.25) 


By applying the chain rule, the spatial derivative of f is described as 
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n—1/2 


n—1/2 
(figura = [(8F /0u)ju lain (2.26) 
The temporal derivative of u can be obtained using Eq. (2.12), 
dpa == (2.27) 
By applying the chain rule again, the temporal derivative of f is derived as 
n—1/2 n—1/2 
(as = [(af/ duu vin (2.28) 


Accordingly, Uz, Fz, Ur, and Fr can be explicitly calculated using the above 
formulas. 

The first-order Taylor expansion in (SEY; is used again to relate Uz“ and Up’ to 
unknowns ui and (ux); at time level n: 


Ax 
U = u’; — z Uo) (2.29) 


* n Ax n 
Ur =u; + Mi (2.30) 


Substituting Eqs. (2.29) and (2.30) into (2.19) and (2.20) yields 


n Ax n At 

Wj — -7 0N} = Ur + (Fr — Fo (2.31) 
n Ax P At 

uj + -g Cj = Urt (Fe — Fr) oy (2.32) 


Adding Eqs. (2.31)-(2.32), one can derive 


1 At 
“= —(U, U — (F -F 2.33 
uj 5 L+ R) + za E R) (2.33) 
and subtracting Eqs. (2.31) from (2.32), one has 


= y= L )+ r ) (2.34) 
— (ux); = U U F F F A 
4 j 2 R L Ax C L R 


Up to now, Eq. (2.33) is presented as an explicit time-marching formula for W’, 
but (ux); still cannot be directly calculated by Eq. (2.34). The reason is because Fc, 
which presents the average value of flux f through the interface CD, has not been 
addressed. In the original CESE scheme, Fc is provided by the Taylor expansion in 
(SE)", i.e., 


2.3 Non-dissipative Core Scheme: a Scheme 17 


A 
Fo= fi- 0H (2.35) 


With a similar procedure to Eqs. (2.25)-(2.28), f ; and ( fj can be expressed 
in terms of u” and (ux). Since uw” is determined by Eq. (2.33), a time-marching 
formula for the only unknown in Eq. (2.34), (ux); , can eventually be derived. 

Note that a complete CESE time step consists of two half steps: a half step 
marching from nodes to centres and another half step marching from centres to 
nodes. The marching schemes for both half steps are identical, except for the indexes 
of the solution points. Usually, the initial conditions of u and ux are specified at the 
initial time, to = 0 and the corresponding boundary conditions for u and uy need to 
be implemented at the boundaries of x. 

The above CESE scheme is referred to as the a scheme [1] in the literature. Define 
a solution vector 


us 
qj=| Ax. | (2.36) 
co (ux); 
and let function f (u) be a linear one: 
f = au, a is a constant (2.37) 


The a scheme, which is mainly expressed by Eqs. (2.33) and (2.34), can then be 
rewritten in a matrix form: 
n n—1/2 n—1/2 
a = Qi) i + Qrad’ i (2.38) 


where the coefficient matrices are 


I[i+v 1—v? 1fi—v -14v? 
a= 3| “4 try f e= » (2.39) 


and v is aconstant, which is defined as 
v=a— (2.40) 


Analogously, we apply Eq. (2.38) to the half-step marching from time level n—1 
to time level n—1/2, and obtain 


qi). = QLA? + Qrati, (2.41) 


n—1/2 n— n= 
drin = Qq; l+ Qraii (2.42) 
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After substituting Eqs. (2.41) and (2.42) into (2.38), a complete time step in the 
a scheme can be formulated as. 


q? = (Q1)’ q7 + QrQr + QeQr)ai! + (Qr)? aiT. (2.43) 


In fact, the numerical dissipation, dispersion, and stability of different CESE 
schemes can be conveniently analysed using their matrix forms like Eq. (2.43). 

An unique feature of this a scheme is its use of Taylor expansion in the inverse time 
direction, to relate the interface flux Fc to the marching variables (see Eq. (2.35)). 
Such a treatment makes the a scheme the space-time inversion invariant and renders 
this scheme non-dissipative [2]. From a historical perspective, the a scheme forms 
the non-dissipative core of subsequent CESE-family schemes. 

The a scheme described in this chapter is restricted to 1D scalar problems and 
uniform meshes. First-order Taylor expansions are used in each solution element in 
the present scheme. Since the a scheme is reversible in time, it cannot be used for prac- 
tical applications, considering that real physical processes are irreversible in time to 
respect the second law of thermodynamics. Accordingly, substantial improvements 
have been made to develop this a scheme into a viable numerical method. In Chap. 3, 
we will describe different approaches to introducing necessary numerical dissipa- 
tion into the CESE scheme, which result in two categories of CESE schemes: the 
central CESE schemes and the upwind CESE schemes. In Chap. 4, we will describe 
the extensions of CESE schemes to multi-dimensional Cartesian and unstructured 
computational meshes. In Chap. 5, high-order Taylor expansions in solution elements 
will be used to achieve high-order accuracy in both space and time, without enlarging 
the spatial stencil or adding time-integration stages. 
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Chapter 3 A) 
CESE Schemes with Numerical E 
Dissipation 


As depicted in Chap. 2, the interface between the two sub-CEs (CD in Fig. 2.7), 
belongs to the SE of (j, n). The flux Fc needs to be calculated through the Taylor 
expansion at point (j, n) toward the inverse time direction. As a result, the a scheme is 
reversible. This violates the second law of thermodynamics. Thus, the non-dissipative 
core suffers from the unphysical oscillations for practical applications. This chapter 
will present the improvements based on the non-dissipative core by introducing 
necessary numerical dissipation. The resulting dissipative CESE schemes can be 
categorized into the central CESE schemes and the upwind CESE schemes. It is 
convenient to begin with the most widely spread CESE scheme, the a—a scheme. 


3.1 a-æ Scheme 


Similar to the a scheme described in Chap. 2, the a—a scheme is also a central scheme 
because of the fact that techniques to calculate upwind numerical fluxes are not used. 
Again, we consider the 1D scalar conservation law of the form 


ðu af(u) _ 


art Ox 


0 (3.1) 


To solve Eq. (3.1) with the a—a scheme, the discretization and the definitions of 
CE and SE keep the same as those for the a scheme, which can be found in Sect. 2.2. 
Basically, the solution algorithm follows that of the a scheme (see Sect. 2.3), except 
for Eq. (2.35) for the flux Fc. Itis worth noting that abandoning the calculation of Fc 
does not affect the correctness of Eq. (2.33) for updating u’. Actually, Eq. (2.33) is 
purely a result of conservation in the entire conservation element (CE)};, irrespective 
of the flux Fc. Therefore, the formula for updating the node value of the unknown u 
remains unchanged, which is written as 
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u" = lu, FUR A er, — FR) (3.2) 

The calculation of Uz, Fz, Ur, and Fr are explicit and has been presented in 
Sect. 2.3, and then uw; can be evaluated using Eq. (3.2). 

The remaining part of the a—a scheme is to update the spatial derivative (ux); 
in a different way from the procedure of the a scheme. It is in this step where the 
necessary numerical dissipation is added into the a—a scheme. First, two different 
estimations for (ux); can be obtained: 


—1/2 1/2 
u" — Bae + (AED) a] 


t); Ax/2 i (3.3) 
ann [Ena + ADU] -u 
(uy dj = Ax/2 (3.4) 


In these equations, the temporal derivatives u; at time level n—1/2 can be readily 
obtained using Eqs. (2.27) and (2.26). Next, (ux); is taken as a weighted average of 
(uy), and (ut): 


(ux), = W (U, UD, a) (3.5) 


where W is a weighted average function with an adjustable parameter a (œ > 0, the 
commonly used values are a = 0, 1 and 2), expressed as 


a —]& 
[xt x + |x xt 


Wera = TS el 


, (3.6) 

The role of this weighted average function is similar to the slope limiter func- 
tions in upwind schemes, and it proves to be effective in suppressing the spurious 
oscillations near a discontinuity. 

The above a—a scheme has a clear and simple logic. However, the dissipation 
of the a—a scheme increases dramatically as the Courant—Friedrichs—Lewy (CFL) 
number (the parameter defined by Eq. (2.40), v = aAt/Ax) approaches zero. In 
practice, for small values of v, the discontinuities in the solution can be smeared 
out. Therefore, the accuracy of the a—a scheme is considered to be CFL-number 
sensitive. 

To demonstrate this shortcoming, we consider f = au (a is a constant) in Eq. (3.1), 
and take œ = 0 in Eq. (3.6). In this case, the weighted average function is reduced to 
a simple arithmetic averaging, and Eq. (3.5) can be written as 


n—1/2 n—1/2 n—1/2 n—1/2 
[wire + A/D a] — [wii + A/D] 
Ax 


us) = » GD 
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which means a central difference approximation at time level n is used to update the 
spatial derivative (ux)j. Recall that in Sect. 2.3 we have shown that a CESE scheme 
applied to the linear scalar convection equation can be written in a matrix form 


qt = Qai Ti + Qrati (3.8) 
for a half step or 
Gj = (Q1)°a} -i + (QLQr + QeQu) aj + Quy aie (3.9) 


for a complete step. Here, the solution vector q is defined by Eq. (2.36) and the 
matrices Qz and Qp are functions of the CFL number v = aAt/Ax. By combining 
Eqs. (3.2) and (3.7), the a—a scheme with a = 0 can also be cast into the form of 
Eq. (3.9) with Qz and Qp as follows: 


Ifit+v 1—v? If i1—v -1+v* 
a=5| he v os] 1/2 -v | po 


Now, consider a limiting case of At = 0 (but Ax > 0), which leads to v = 0. It 
can be shown that Eq. (3.9) leads to 


1f 1/2 1 a 1/30] ,-,, 1|1/⁄2 -1 i 
n— L n Z n = n .11 
a T aleia aC Fala 12/41 61D 


However, a reasonable time-marching scheme should guarantee that q; = qi! 
when At vanishes, which indicates that the a—w scheme suffers from numerical error 
when the CFL number is very small. Further numerical experiments show that, in 
practice, when the CFL number v becomes smaller than 0.1, the excessive numerical 
dissipation of the a—a~ scheme can be remarkable. 


3.2 Courant-—Number-Insensitive Scheme 


To overcome the shortcoming of the a—a scheme, a Courant-Number-Insensitive 
(CNI) scheme was constructed [1]. This improvement is based on two requirements: 
(1) the CNI scheme should reduce to the non-dissipative a scheme when the CFL 
number v = 0 and (2) the CNI scheme should resemble the a- scheme when 
|v| = 1. In this section, the construction of the CNI scheme will be present, and then 
the relationship between the CNI scheme, the a scheme, and the a—w scheme will be 
shown. 

First, the formula for updating the node value of the unknown u remains the same 
as Eq. (3.2), which has been used in a scheme and the a—a scheme. By substituting 
formulas for Uz, Fz, Ur, and Fr (Eqs. (2.21)—(2.28)) into Eq. (3.2), the explicit time 
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marching formula for u’; can be written as 


n 1+v yr 1— n-1/72, 1-V n1/2 n-1/2 
uj = ete "Axtu Mx) jp t = a 3 T axtudintf 
(3.12) 


Then, an algorithm for updating (ux); is devised. As shown in Fig. 3.1, the points 
(j — 1/2, n) and (j + 1/2, n) are denoted by V_ and V4, respectively. The point M_ 
is the midpoint between (j — 1/2, n) and (j, n), and M, is the midpoint between (j 
+ 1/2, n) and (j, n). We define two important points P_ and P+, which are on the 
line segments V_M_ and V,M,, respectively. The distances between P and node 
(j, n) are marked in Fig. 3.1. It is clear that P+ coincides with V+ when v = 1, and 
they approach Mz as v — 0. The values of u(x, t) at P4 are estimated by Taylor 
expansion at points (j + 1/2, n — 1/2): 


no A (1 — |v) Ax ü= 
u(Ps) = uzi + + Zu Da F —q — Haji (3.13) 


where igs and (u x) js as are known, while Ca as can be obtained using Eqs. 


(2.27) and (2.26). 
Next, u(P_) and u(P,,) are used to construct two different approximations of (ux); 


ui —u(P_) 


(a 3.14 
(Gi); = TFAA ii 
and 
u(P,) — u” 
at)" = ———__/_ 3.15 
6); = TFAA ai 
V P. M. (J, n) M, P, V, 
- n 
(I-\v|)Ax  (Hv])Ax (1+|v|)Ax 
= 4 
2 
l, 
x 
as l i | l 
[i-3. n- J Lie. n-5] 


Fig. 3.1 Definition of points in CNI scheme 
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Finally, the CNI scheme use a weighted average of (ay; and Qty; as the updated 
value of (ux)}, in a way similar to the a—a scheme. However, the weighted average 
function in the a—a scheme, i.e. Eq. (3.6), is replaced with a more sophisticated one 
as 


[1+ fdvdo jar + [1+ rade ]@oy wie 
ux . l l 
j 2+ Fv)[ 619 + HF] 


where 
Cra 
(s+); = - - (3.17) 
min(|@y"], |@24)) 
and 
flv) =0.5/ 1v]. (3.18) 


Hence, the overshoot phenomenon near discontinuities can be suppressed by the 
artificial dissipation, just like the a—w scheme. Moreover, the dissipation brought by 
Eq. (3.16) can be adjusted dynamically according to the CFL number v. 

When v = 0, we can show that the CNI scheme will reduce to the non-dissipative 
a scheme. Substituting v = 0 (which also means At = 0) into Eqs. (3.12) and (3.13) 
yields 


1 a12 1 -1/2 l aap 1 -1/2 
ui = z“ + z AU) + z“ j+1/2 — g ATU) j (3.19) 


and 


n— Ax TE 
UCP) = uipa F g U a (3.20) 
From Eqs. (3.19) and (3.20), we can get 


g ei ap 


By using v = 0 and Eq. (3.21), it is readily shown that 


n P,) —u(P_ a)n 
ED = E = a (3.22) 


Thus, the weighted average of (û7 yj, and (ay) is equal to [u(P4) — u(P_)|/(Ax/2), 
which gives 
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n 1 n—1/2 Ax n—1/2 n—1/2 Ax n—1/2 
ux); = |e Gq U = Uj. t g Uj- (3.23) 
Recall the expression of the a scheme (see Eqs. (2.38) and (2.39) in Chap. 2) and 

let v = 0, we can reproduce Eqs. (3.19) and (3.23). Therefore, in the limiting case 

of v = 0, the CNI scheme and the a scheme are identical. 
As for the other limiting case Ivl = 1, Eqs. (3.13)-(3.15) lead to 


a MEH [wena + AHO AT 
(az)" = (3.24) 
J Ax /2 


and 


n—1/2 n—1/2 n 
n [wine + ADU | -e 
(ât) = L HeJ j (3.25) 
J Ax/2 


There is no difference between AE) in the CNI scheme and (u>)" in the a—a 


scheme. Hence, the algorithms for updating (ux); in the CNI and the a—a schemes 
are basically the same when lvl = 1, except for the specific forms of the weighted 
average function. 

Both the a—a scheme in Sect. 3.1 and the CNI scheme in this section belong to 
the central CESE schemes. Extensions of these schemes to 2-D and 3-D cases for 
various systems of conservation equations (e.g. Euler equations for compressible gas 


dynamics) are straightforward and have been well implemented. 


3.3. Upwind CESE Scheme 


The aforementioned central CESE schemes are used to solve nonlinear hyperbolic 
systems of conservation laws. However, they did not explicitly resort to the knowl- 
edge of characteristics or eigenvalues of the systems. As an alternative approach, 
the characteristic-based upwind CESE scheme proposed by Shen et al. [2], elegantly 
combines the basic ideas of the CESE method and the upwind numerical flux tech- 
nique in the Godunov-type FVM method. The upwind CESE scheme is naturally 
CFL-number-insensitive, which means it does not suffer from the drawback of the 
a-a scheme presented in Sect. 3.1. For some challenging CFD problems such as the 
simulations of detonations and multiphase flows, the upwind CESE scheme captures 
discontinuities in flow fields with improved accuracy and robustness, especially for 
contact discontinuities (e.g., material interfaces). 
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3.3.1 Construction of Upwind CESE Scheme 


For illustrative purposes, the 1D scalar conservation law (Eq. (3.1)) is considered 
again. Before introducing the upwind CESE algorithm, some elementary concepts 
about the space-time discretization need to be clarified. The same computational 
mesh and solution points as shown in Figs. 2.3 and 2.4 are adopted. The definition of 
CEs is also retained (see Fig. 2.5). Nevertheless, some modification is made to define 
the SEs, as sketched in Fig. 3.2. In comparison to the (SE) in Fig. 2.6, the newly 
designed (SEY; no longer contains any part that allows for t < t,. Apparently, the SEs 
pave the whole space-time domain without overlap. It is assuemed that u(x, t) and 
f (, t) are piecewise linear. For instance, inside (SE), they can be approximated 
by the first-order Taylor expansion at point (j, n). The boundaries of conservation 
element (CE) belong to three solution elements: DE and DF belong to (SE ts AC and 


AE belong to (SE) ips and BC and BF belong to GE. The interface between 


(GE) 12 and Gas i.e. line segment CD, splits the conservation element (CE)’ 
into two sub-CEs: (CE); and (CE*)¥. 

The construction of the upwind CESE scheme closely follows the framework 
of the non-dissipative core scheme (i.e., the a scheme in Sect. 2.3). Indeed, Eqs. 
(2.16)-(2.34) still hold for the upwind CESE scheme, and they can be derived by 
the same procedures as presented in Sect. 2.3. Therefore, just like the a scheme, the 
a—-a scheme, and the CNI scheme, the formula to update u’; in the upwind CESE 
scheme remains the same as Eq. (3.2). However, for the purpose of updating (u,)’,, 
the Eq. (2.34) will be utilized, which is the direct result of the space-time integral 
form of conservation law. For convenience, here we recall this equation: 


es tie Ore ae) (3.26) 
4 Mi = 5 R L FAX C L R . 


Note that this equation is actually discarded in the a—a and the CNI schemes, but 
both the a scheme and the upwind CESE scheme make full use of it. Furthermore, 


1 l 
Uan) 


Fig. 3.2 Definitions of CE and SE for upwind CESE schemes 
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the key difference between the a scheme and the upwind CESE scheme lies in the 
treatment of Fc in Eq. (3.26), which denotes the average flux through the interface 
CD between (CE); and (CE+); (see Fig. 3.2). 

As can be inferred from the new definition of solution elements in Fig. 3.2, the 
operation used in the a scheme to link Fç with (ux yi is forbidden in the upwind CESE 
scheme. Instead, we intend to obtain Fc from the known data at time level n — 1/2, 
and then insert Fc into Eq. (3.26) to get (ux);. Toward this end, a local Riemann 
problem can be built at the midpoint of CD (denoted by (j, n — 1/4)), because CD is 


the interface between (SE); i and (SE) 712 To be specific, by Taylor expansion 


in (SE) i the left state at point (j, n — 1/4) is evaluated as 


n—1/4 Ax n 2 At n 2 
(uL); = U, + F Uo + g UD (3.27) 
n—1/2 


Meanwhile, Taylor expansion in (SE);,\/5 provides the right state at point 
(G, n — 1/4): 


ma Ax ñ At ye 
(ua) = Ur = T U + y UV (3.28) 


Recall that Uz and Up are given by Eqs. (2.21) and (2.23). According to the defi- 
P ‘ n—1/4 n—1/4 . 
nition of solution elements, the values of (uz) j and (ur) j are not necessarily 


equal to each other, and in general (u D” * and (u a 4 form the discontinuous 
initial data on the left and right sides of the “diaphragm” CD. 

Next, the average flux Fc through CD can be evaluated based on the local Riemann 
problem with initial data (3.27) and (3.28) as 


Fo = P( (ui), (ue) ") (3.29) 


where F stands for any upwind numerical flux solver, or loosely referred to as 
Riemann solver. Consequently, the special derivative (ux); can be updated explicitly 
by Eq. (3.26). It is worth noting that, in the upwind CESE scheme, the time marching 
formula Eq. (3.2) for u’; has no concern with the upwind procedure above, since flux 
Fc never appears in Eq. (3.2). 

In case of strong discontinuities, using proper limiters for the derivatives in Eqs. 
(3.27) and (3.28) bocomes crucial to suppress spurious oscilations. When limiters 


are used, the reconstruction of (u D” and (u R" * is written as 
= Ax At 
Up = UL + Su + uy, (3.30) 
z A At 
ie = Ur — >u? + u? (3.31) 
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where ut, u®, uř, and uř are the limited derivatives. In this book, the weighted 


biased averaging procedure limiter (WBAP-L2) [3] is adopted. The limited slopes 
(spatial derivatives) are 


a-1/2 
= (ux) wd, 0}, 0F), OF = (Ur Ua NEE) oe (ux) i2 


j-1/2 n—1/2 I2 = n—1/2? 
(Ux) j-1/2 (Ux) j-1/2 
(3.32) 
sein E —U1)/(Ax/2) pp Us) 
ue = (Ux) WC, OF, 07), OF = naif 82 = ap? 
(Ux) j41/2 (Ux) 41/2 
(3.33) 
where the limiter function is 
S+V/O+/ ifo andô, > 0 
WA, 61,05) = 4 SHIFTER OE (3.34) 
0, else 


Once the limited slopes are obtained, the limited temporal derivatives can be 
calculated by the chain rule and the conservation equation itself, i.e., 


ub = L of 


tee (3.35) 


Uy 
u=UL ðu u=UR 


t -u a 


3.3.2 Scheme for Linear Scalar Convection Equation 


Let f = au, then Eq. (3.1) becomes the linear scalar convection equation. Without 
loss of generality, a is assumed to be a positive constant. For this simple situation, the 
exact solution of the local Riemann problem with initial data (u D” 4 and (u R); HA 
can be readily obtained. If the limiter is not used, the interface flux Fc (Eq. (3.29)) 
is 


Fe salu +S (u); A we At (a > 0) (3.36) 


With this Fc, the upwind CESE scheme, which is a combination of Eqs. (3.2) 
and (3.26), can be written in a matrix form as Eqs. (3.8) and (3.9), in which g = 
[u, (Ax/4)u,]’ and the matrices Qz and Qpr are functions of the CFL number v = 
aAtlAx: 


1] 1+v 1-—v? I[1-v -1+v? 
= — = — .37 
Q | (ap -1+4 aK Tien (an 
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n—1/4 n—1/4 


If the initial data (uz) and (up) are reconstructed using the limiter, the 
matrix form shown as Eqs. (3.8) and (3.9) still holds, but the matrices Q; and Qpr 
should be 


If l+v 1—v? 1fi—v -14v? 
Q=- 2 h Qr=5 2 
2| —1 +v -14+20+4+ L)» + (1 — 2¢z)v 2L1l—v -1 +v 
(3.38) 
where 
(Ur — Uz)/(Ax/2) Gas 
u= WU Of, 0f),0f = 4 oy = 453.39) 
(Ux) j-1/2 (Ux) j-1/2 


The scheme without limiter Eq. (3.37) can be regarded as the special case of 
Eq. (3.38) with the limiter function ġz = 1. 

Based on matrices Qz and Qp in Eq. (3.38), we can examine the property of the 
upwind CESE scheme at the limiting case of v = 0. Recall that the a—æ scheme 
becomes very diffusive when v — 0 and cannot guarantee that qi > q; as 
At — 0 (but Ax is a finite constant). However, the upwind CESE scheme does not 
suffer from such a deficiency. A direct evaluation of matrices Qz and Qp in Eq. (3.38) 
as v — 0 shows 


2 2 
if 1 1 Ifi —1 
Zaga = 2 ay — 
Qz TA 1O se i 


Ii 1 1 -l 1/1 -1 1 1 
QLQr + QRQ: > TE zalk HERI mle “| = 


Consequently, the time marching scheme for a complete time step (Eq. (3.9)) can 
always ensure that q —> qi! as At — 0, regardless of the form of the limiter 
function. Unlike the a—a~ scheme, the upwind CESE scheme is inherently CFL- 
number insensitive. 


(3.40) 


3.3.3 Scheme for Euler Equations 


The 1D unsteady Euler equations for a perfect gas are written as 


dU OF 
a Al 
ot Ox i ie 


where U = [p, pu, EJ’, F = [pu, pu? + p, (E+ pu] are the vectors of the 
conserved variables and the inviscid flux. Here, p, u, and p represent the density, 
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velocity, and pressure, respectively. The total energy per unit volume is denoted by 
E, which is E = p/(y — 1) + pu*/2 with constant ratio of specific heat y = 1.4. 

The upwind CESE scheme for a scalar equation can be directly applied to solve 
each component of U. However, the upwind procedure to calculate the flux through 
the interface between (CE )’ and (CE*); needs a substantial extension, because of 
the complexity of the eigen-structure of Eq. (3.41). Still, the flux vector Fc can be 
expressed as 


Fe = Ê( UD, UD), (3.42) 


where F stands for any appropriate upwind numerical flux solver, e.g., approximate 
Rimann solvers and other upwind flux functions. Moreover, (Uz i. 1/4 and (U R)j 1/4 
are the data on the two sides of line segment CD in Fig. 3.2, which are reconstructed 
with an appropriate limiter and serve as the discontinuous initial data of the local 
Riemann problem. Fortunately, any efficient and robust upwind flux solver [4, 5] 


exsiting in the FVM literature can be employed to calculate Fc. 


3.3.4 Remarks on Upwind CESE Method 


For the purpose of capturing discontinuities, the upwind CESE method introduces 
necessary numerical dissipation by the upwind procedure to calculate Fc (see Eqs. 
(3.29)-(3.35)). This fact might lead to confusion between the classical upwind FVM 
[4, 5] and the present upwind CESE method. In the former, the upwind flux technique 
to tackle local Riemann problems is viewed as the building block of the whole method. 
In the latter, the upwind flux technique plays a different role. 

To shed light on this issue, the interval [j — 1/2, j + 1/2] is taken as a repre- 
sentative control volume to analyse, which is called “cell” in the FVM, and marked 
as line segment AB in Fig. 3.2. In this section, three fluxes are related with [j — 
1/2, j + 1/2], denoted by Fz, Fr, and Fc. Apparently, these fluxes can be classi- 
fied as the flux through the boundary of the control volume (Fz and Fp) and the 
flux through the diaphragm inside the control volume (Fc). For the upwind CESE 
scheme, the upwind procedure is performed only to calculate Fc, while fluxes through 
the boundaries are not linked to any local Riemann problems. Owing to the time- 
marching strategy on the staggered mesh, information at points x; 1/2 and xj41/2 are 
ready before evaluating Fz and Fr. Thus, the procedure to calculate Fz and Fg only 
involves Taylor expansion inside the SE, the definition of physical flux function f (u), 
and the Cauchy-Kowalewski procedure (steps to derive u; and f, from ux as shown 
in Eqs. (2.26)—-(2.28)). On the contrary, the non-staggered FVM employs the upwind 
flux solver to get all fluxes through boundaries of control volumes, because each Fz 
or Fr should be treated as the result of a local Riemann problem. Usually, Fc is not 
considered in the FVM. 


32 3 CESE Schemes with Numerical Dissipation 


Recall that Fc is absent in Eq. (3.2). Therefore, the upwind procedure in the CESE 
scheme never affects the time-marching algorithm for the average value of u(x, t) on 
li — 1/2, 7 + 1/2], but solely redistributes u(x, t) on [j — 1/2, j + 1/2] through the 
calculation of u, using Eq. (3.26). However, the upwind FVM relies on two upwind 
numerical fluxes fs 1/2 and fa 1/2 to update the average value of u(x, t) on [j — 1/2, 
j + 1/2]. This is the reason why we claim that the upwind CESE method utilizes the 
upwind technique in a different way from the traditional upwind FVM. 


3.4 Comparison of Different CESE Schemes 


So far, four CESE schemes, namely the a scheme, the a—a scheme, the CNI scheme, 
and the upwind CESE scheme, have been introduced. To make a comparison, their 
main formulas and properties are listed in the following tables. Table 3.1 indicates that 
the CNI and upwind CESE schemes are more favourable for the practical simulations, 
since they can capture discontinuities with reasonable numerical dissipation and 
they are free of the sensitivity issue which is encountered by the a—a scheme. The 
key formulas in each scheme are shown in Table 3.2. It is notable that all CESE 
schemes share a common approach to updating u’. Only the a scheme and the 
upwind CESE scheme update (ux yj are based on the conservation law for sub-CEs. In 
contrast, the a—a and the CNI schemes construct (ux); using the finite-difference-like 
approximation followed by some kind of weighted averaging technique. 


Table 3.1 Properties of four different CESE schemes 


Central or upwind Dissipative or CFL-number sensitive or 
non-dissipative insensitive 
a scheme Central Non-dissipative Insensitive 
a—a scheme Central Dissipative Sensitive 
CNI scheme Central Dissipative Insensitive 
Upwind CESE Upwind Dissipative Insensitive 


Table 3.2 The formulas to update ui and (ux); in each CESE scheme 


Formula for u’ 


Formula for (ux)j 


a scheme 


a—a scheme 


CNI scheme 


a 
uj = 5 UL + UR) 


+^ (r Fr) 
2Ax z s 


Upwind CESE 


£2 (ux) = (Ur — UL) + Fc — Fi — Fr) 
u= whas WEY, a) 
Ui = WAGO GY, v) 
EF ux)" = 4 (Ur — UL) + (Fe — Fi — Fr) 
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Table 3.3 The matrices Qz and Qp in Eq. (3.8), which is the matrix form of the CESE method 
applied to the linear scalar convection equation (Eq. (3.1) with f = au, a > 0) 


a scheme z 


2 2 
a—a scheme («œ = 0) | fe] 


Upwind (without limiter) Z 


l+v 1-—v? | 


When applied to the linear scalar convection equation (Eq. (3.1) with f = au and 
a is a positive constant), it is possible to rewrite the CESE scheme in a very compact 
matrix form. In Table 3.3, the coefficient matrices for different CESE schemes are 
listed. The matrices for the CNI scheme are not tabulated due to their lengthiness. 
Such Qz and Qpr, together with Eqs. (3.8) and (3.9) are useful for revealing the 
instrinct properties of different CESE schemes, e.g., the stability of each scheme. 


3.5 Numerical Examples 


Here, the 1D scalar convection problem and Sod’s shock-tube problem with uniform 
grid size of 0.01 are computed using different CESE schemes. The corresponding 
C + + source codes implementing the a—a CESE scheme are presented in the 
Appendix. For Eq. (3.1) with f = u, a square wave propagates in a computational 
domain [—1, 1] till t = 2.0 with the initial condition described as 


1, if -0.5<x<0.5 
CD ea jse (3.43) 


The periodic boundary condition is imposed on both ends. The results are plotted 
in Fig. 3.3. For CFL = 1, all three schemes provide a virtually exact result. When 
CFL = 0.8, both the a-a and the CNI schemes suffer from strong oscillations on 
the upwind sides of the discontinuities, while the upwind scheme provides a result 
with satisfactory accuracy. If CFL is further decreased to an extremely small value, 
the CNI and upwind scheme can capture the discontinuities very well, but the large 
dissipation in the a—a scheme abnormally smears out the discontinuities. This case 
proves the CFL insensitivity of the CNI and upwind CESE schemes. 

In the Sod’s shock-tube problem [6], the gas is initially separated at x = | with 
left and right states (p, u, p) = (1.0, 0.0, 1.0) and (p, u, p)r = (0.125, 0.0, 0.1). 
In Fig. 3.4, the computed density profiles at t = 0.4 are shown. For CFL = 0.8, all 
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(a) a-a CESE (b) CNI CESE 


—cFL=1 
— CFL=0.8 
—+— CFL=1E-4 


-1.0 OS 0.0 0.5 10 
x 


(c) Upwind CESE 
Fig. 3.3 Solutions for scalar convection equation: square wave problem at t = 2 
the schemes can capture the wave structures precisely. For an extremely small CFL 


number, the CNI and upwind schemes maintain great accuracy, but the a-æ scheme 
seriously smeared out the rarefaction, contact, and shock waves. 
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Fig. 3.4 Density profile of Sod’s shock-tube problem at t = 0.4 


References 


1. Chang, S.C, & Wang, X.Y. (2002). Courant number insensitive CE/SE Euler scheme. In 38th 
AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit. 

2. Shen, H., Wen, C. Y., & Zhang, D. L. (2015). A characteristic space-time conservation element 
and solution element method for conservation laws. Journal of Computational Physics, 288, 
101-118. 

3. Li, W., & Ren, Y. X. (2012). The multi-dimensional limiters for solving hyperbolic conservation 
laws on unstructured grids. Journal of Computational Physics, 230, 7775-7795. 

4. Toro, E. F. (2009). Riemann solvers and numerical methods for fluid dynamics: A practical 
introduction. Springer. 

5. Kitamura, K. (2020). Advancement of shock capturing computational fluid dynamics methods: 
Numerical flux functions in finite volume method. Springer Nature. 

6. Sod, G.A. (1978). A survey of several finite difference methods for systems of nonlinear 
hyperbolic conservation laws. 


36 3 CESE Schemes with Numerical Dissipation 


Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 
International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, 
adaptation, distribution and reproduction in any medium or format, as long as you give appropriate 
credit to the original author(s) and the source, provide a link to the Creative Commons license and 
indicate if changes were made. 

The images or other third party material in this chapter are included in the chapter’s Creative 
Commons license, unless indicated otherwise in a credit line to the material. If material is not 
included in the chapter’s Creative Commons license and your intended use is not permitted by 
statutory regulation or exceeds the permitted use, you will need to obtain permission directly from 
the copyright holder. 


Chapter 4 ®) 
Multi-dimensional CESE Schemes gaa 


The previous chapter has shown that the necessary numerical dissipation can be 
introduced in 1D CESE schemes through either a central or upwind approach. 
In this chapter, we present the extensions of 2D CESE schemes based on carte- 
sian meshes, followed by a detailed description of implementation on unstructured 
meshes. Both the central schemes and upwind schemes will be presented. Then, 
numerical examples will be provided to demonstrate the capabilities of the present 
schemes. 


4.1 CESE Schemes on Cartesian Meshes 


4.1.1 The Improved a- CESE Scheme 


The marching scheme requires physical variables and their spatial derivatives at each 
mesh point. In the pioneer development of the 2D CESE solver [1], the domain is 
discretized by congruent triangles. Based on a similar technique, Zhang, Yu, and 
Chang [2] reported a further extension of the 2D CESE scheme on quadrilateral 
meshes. In the above versions of CESE schemes, the solution points are solely 
updated at the cell centers with a staggered stencil. As a result, a generalized flux 
technique was also proposed as a post-marching procedure to handle the “flux decou- 
pling” problem when two mesh points cohost one sub-CE. Alternatively, an improved 
2D CESE scheme [3, 4] was proposed with a new definition of SE and CE. In this 
updated scheme, the entire space-time region is divided into non-overlapping CEs, 
making it convenient and simpler for calculation and straightforward for extension 
to 3D scheme. Figure 4.1 shows the space-time geometrical configuration of the 
improved 2D CESE scheme with a uniform rectangular mesh. The spatial projec- 
tions of the solution points at (n —1/2)th and nth time levels are denoted as W and 
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@, respectively. The solution is updated alternatively between the cell centers @ and 
cell vertices W. Here At/2 = t, — t,-1/2 and P, P’, and P” subsequently repre- 
sent the points at three successive timesteps. The conservation element of point P’, 
denoted by CE(P’), is defined by the space-time hexahedron ABCDA’B’C’D’, i.e., 
the union of four sub-CEs: CE”? (P’), CEŁ" (P°), CE®?(P’), and CEFY (P’). In partic- 
ular, CE? (P>) is defined as AFPEA’F’P’E’. Other sub-CEs are defined similarly. 
Next, consider the solution element of P’, SE(P’), which consists of three orthogonal 
planes, A’ B’ C’ D’, EGG” E”, and FHH” F”. The physical flux vector is assumed to 
be smooth within each SE and can be approximated by Taylor expansions about the 
mesh point associated with the SE. In the following, we present the construction of 
the 2D CESE solver using similar techniques as described in the 1D scheme. 


(c) SE(P') (d) Fluxes in the sub-CEs 


Fig. 4.1 Grid points in the spatial domain, definitions of CE and SE, and corresponding fluxes in 
a CE for the 2D CESE scheme on rectangular meshes 
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In this section, we consider the 2D scalar hyperbolic conservation law 


au of ə 
ee! ee (4.1) 
ot ox ody 


By using the Gauss’ divergence theorem, it is shown that Eq. (4.1) can be written 
in the form of 


$ h-ndS=0, (4.2) 
S(CE(P’)) 


where h = (f, g,u) is the space-time flux vector, S (CE(P’ )) is the boundary of 
CE(P’), n is the unit outward normal vector on the surface of the control volume. 
Note that the u, f, and g are approximated by first-order Taylor expansion about point 
P’. For any (x, y, t) € SE(P’), let 


u(x, y,t) = u(x, dy, dt) p, (4.3) 
f(x, y,t) = f(x, dy, dt) p, (4.4) 
g(x, y, t) = g(dx, dy, dt) p, (4.5) 


where X (ôx, dy, dt) y denotes the first-order Taylor expansion about point N as 
X (dx, dy, dt)y = Xy + (Xx) Ny bx + (Xy) yoy + (X) vot, (4.6) 
and 
ôx =x — xy, ôy = y — yy, ôt =t — ty. (4.7) 
In addition, the spatial and temporal derivatives of f (u) and g(u) can be derived 
by the chain rule (Eqs. (2.26) and (2.28)) and u; = — fy — gy. 


Meanwhile, one obtains local flux conservation of the four sub-CEs by integrating 
Eq. (4.2) over their surfaces, given by 


, AxAy Ax Ay AyAt Ax At 
Urn = Urp + (Frp — Fcp) + (Gpr — Ger) , 
4 4 4 4 
(4.8) 
, AxAy Ax Ay AyAt Ax At 
Ut r- Urp 4 * (Fcp — Fro) ar (Gpr — Ger) a. 


(4.9) 
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, AxAy Ax Ay an AZAT 
Uku = = Ur + (Feu — Fru) + (Ger — Gur) 
m 10) 
, AxAy Ax Ay ala 
Uly = Uzy + (Ftu — ppe an (Gcr — Gut) 
4 4 
ar 11) 


where Uzp, U’:p, Frp, Fcp, GpL, and Gcz are the average fluxes on the surfaces 
of CE (P’), AFPE, A’F’P’E’, AEE’A’, FPP’F’, AFF’A’, and EPP’E’, respectively. 
Similar definitions are applied to the other sub-CEs. By definition, the fluxes U;p, 
Fp, and Gp, can be determined via the Taylor expansion in SE(A). Itis also important 
to emphasize that the fluxes Fcp and Gcz that denote the fluxes across the “inner” 
surface of two neighboring sub-CEs, however, are not trivially available due to the 
presence of discontinuities. In fact, it will be shown later that all the fluxes across the 
interfaces among sub-CEs vanish in the formulation of the time marching scheme of 
u. Apart from this, based on the assumption that the distribution within each SE is 
linear, the average fluxes on the exterior surfaces of CE/?(P’) read 


Ax A 
Uio =u( $E, 2.0) (4.12) 
A 
REEN ON 4.13 
Lp 7“ et p’ (4.13) 
z 
Ax At 
T e 0, — (4.14) 
se aa a 
Ay At 
Kase(0,—.— (4.15) 
4° 4), 


Fluxes are expressed similarly for the other three sub-CEs. 

To proceed, by adding Eqs. (4.8)—(4.11), the fluxes through the interfaces between 
sub-CEs are well balanced. Consequently, the value of the solution point P’ can be 
derived as 


u(P’) = {Uo + Uro + Uru + Uw) + £ = (Fp Frp — Fru + Fru) 


MEA P + Grp — Gru — Gru) 
4Ay 
(4.16) 
Moreover, by replacing the average fluxes in Eq. (4.16) with the Taylor expansions 


about their corresponding SEs centers, we can obtain the explicit time marching 
scheme for u(P’) as 
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u(P’) = la+ f+ F, (4.17) 


where 


a=u 2x EY Hu Ax A o Hu ca AY o +u >A AY 9 3 
4° 4°), r 4° 4 Je 4? a p 
Ay At Ay At Ay At 

FLO, , f{9, =, F{ 9, A , 
4° 4)p 4° 4), 4° 4 je 


Ax At Ax At Ax At Ax At 
8g 0, g .0, 8 »0, g 0, i 
4 4/4 4 4)p 4 4 jc 4 4 jp 


With respect to the spatial derivatives, by using the continuous assumptions at 
points A’, B’, C’, and D’, the following relations can be formulated 


SI 
ae 
F R 
S 
> 
< 
> 
< 
ae 


os| 


Ax Ay At 
u{ -—,-—,0]} =uļ|0,0,— |], (4.18) 
2 2 p 2/4 
Ax Ay At 
u{ —,-—,0) =u{0,0,—}] , (4.19) 
2 2 ; 2 Jpg 
Pp 
Ax Ay At 
u{ —,—,0] =u(0,0,—] , (4.20) 
2 2 p 2 jec 
Ax Ay At 
u{ ——-, —,0} =u {0,0,—] . (4.21) 
2 2 p 2 Jp 


After mathematical manipulation of Eqs. (4.19) and (4.20), (4.18) and (4.21), 
(4.20) and (4.21), (4.18) and (4.19), respectively, the spatial derivatives at solution 
point P’ can be explicitly calculated as 


T(P’) = A fa oo >) + (0 0, ) —2uP’) (4.22) 
va ~ Ax n(o. a) B ae 2 jec “ l 
(P) = ' (0 i ‘) (0 : ‘) APD Ea 
"x =-5 | 2 Ig A ie ~ | 
= l 0,0 ") + (0 0, of) a (P’) (4.24) 
E Ay a(o. ae pu D í ' l 


-Ph =- (0 0 >) (0 0 >) 2u(P’) (4.25) 
uy =-5 hu oan oti re oe | š 
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One note that, two derivatives were estimated along both spatial directions 
(denoted by superscripts + and —). Recall that to a weighted averaged function 
Eq. (3.6) is again adopted to suppress numerical wiggles, such that 


il P= W (Cie) wis (te ) po @); (4.26) 


uy(P’) = W((ws) ps (ut) p» a). (4.27) 


The Eqs. (4.17) and (4.22)-(4.27) constitute the first half-step marching scheme 
from cell vertices to cell centers. The second half-step scheme that marches from 
centers to vertices follows a similar approach. As described above, constructing the 
improved 2D a—a CESE scheme based on rectangular grids is genuinely simple and 
could be straightforwardly extended to the 3D scheme. 


4.1.2 CNICESE Scheme 


When applying the above 2D a—a schemes [1—4] for solving problems with highly 
non-uniform meshes, the local CFL number may vary significantly from its stability 
limit of 1 to approaching 0. Namely, the solutions tend to smear out by decreasing 
the CFL number. The Courant number insensitive (CNI) scheme was proposed to 
address this issue to alleviate the sensitivity problem by improving the strategy in 
updating the derivatives [5]. A new point, Im, is designated dynamically moving 
along the line segment connecting the vertices and the centers M’m of each sub- 
CEs (Fig. 4.2). Analogous to the 1D CNI scheme (Sect. 3.2), if we define the local 
and global Courant numbers as v and vo, the interpolated points are designated to 
approach the centers when v/vg — 0, and approach the cell vertices when v/v9 — 1. 
For example, the coordinates of 7 í are 


xI) = Ha) + (1 = ~) scan), (4.28) 


yl) = —y(A’) + (1 = z). (4.29) 
vo vo 


If we approximate u(I D) at t = tn from the expansions about A at t = tn—1⁄2 and 
P’ att = tn as 


At 
u(x, y+, 0p = u(i) = u(x. ôy, >) (4.30) 
A 


where 5x4 = x(I,) —x(P’), ô&y+ = y(U,) — y(P’), 6x_ = x(I;) —x(A), and dy_ = 
yd D) — y(A). Similar relations can be derived for each Iy: Combining Eq. (4.30) 
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Fig. 4.2 Definition of points 
in 2D CNI CESE scheme 


for u(I i) and that for u (1 a) one can explicitly formulate two independent conditions 
for calculating the spatial derivatives of u(P’). Repeating the same procedure for 
combinations of u(1;) and u(13), u(1;) and u(15), u(1;) and u(14), we arrive at four 
sets of spatial derivatives. The weighted average function is then used to calculate 
the optimal derivatives for capturing the discontinuities. 


4.1.3 Upwind CESE Scheme 


In the 2D upwind CESE scheme that Shen et al. [6] proposed, the CE retains the 
same form as the central CESE scheme. The time marching scheme of u is the 
same as in the a-a~ scheme. The upwind procedure only affects the calculation of 
the spatial derivatives, and global space-time conservation is ensured. In analogue 
to the definition in the 1D upwind scheme, for example, SE(P’) is defined as cuboid 
A’B’C’D’A”B”C”D” (Fig. 4.3), where the physical flux vector is approximated by 
the Taylor expansion about point P’. The average values at the top surfaces of the 
four sub-CEs of CE(P’) read 


Ax Ay 
Vin=u(-F.-2.0) | (4.31) 

F Ax Ay 
Uko =u( $£, -2.0) , (4.32) 


; Ax Ay 
Uru =ul —-,—,0) , (4.33) 
4 ji 
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At 
2 
x 
Ax/2 
Fig. 4.3 Definition of SE of the upwind CESE scheme 
U! ar S (4.34) 
=z ul == ; ; 
Ai 4'4’ Jp 


Referring to Fig. 4.1d, by substituting Eqs. (4.31) in (4.8), (4.32) in (4.9) and 
subtracting the two, one obtains 


uP (P) 2 [veo ULp) + a (2Fcp -Ftp — Frp) + T (Cpr Gcr-GpL+ Gcu) | (4.35) 

Equation (4.35) represents the spatial derivative in the x direction based on the 
information of the lower pair of CEs. Note that in contrast to central schemes, the 
upwind scheme requires additional estimations of inner fluxes (e.g., Fcp, Gcr, and 
Gcı in Eq. (4.35)) to compute spatial derivatives. Analogously, with the aid of Eqs. 
(4.31) and (4.34) and by combining with Eqs. (4.10) and (4.11), Eqs. (4.8) and (4.11), 
Eqs. (4.9) and (4.10), respectively. The following equations for spatial derivatives 


can be formulated 


2 At At 
ul (P) Te [Uru Uru) Aa (2Fcu — Fru — Fru) 4 Dy (Gcr - Gur- GcL +4 Gus)], (4.36) 
«LoP = 2| (Uru Urp) + ~ (Fry -Fey -Frp + Fep) + tee etna Co) |. (4.37) 
y Ay | VLU — VLD) + Ky ELU u — FLD D) + ay OCL- OpL-GuL) f 

2 t At 
uy (P’^) [Uru Urp) Ay (Feu Fru — Fcp + Fro) Ay COCR GDR Gur)}, (4.38) 


Ay 
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The final spatial derivatives are calculated by using WBAP limiter to remove small 
noises near the strong discontinuities, e.g., 


ux (P!) = ul (P'WBAPEZ(1, 01, 0) = [RENY eh] wearta, 01, 62). (4.39) 
where 


6, = uP(P’)/ul (P’), 


6. = ul (P')/uC(P’). 


In the expressions for spatial derivatives in Eqs. (4.35)—-(4.38), the fluxes on CE(P’) 
surfaces can be easily approximated by the Taylor expansions of the values from t 
= ftn—1/2. To complete Eqs. (4.35)-(4.38), the fluxes through the inner boundaries 
of CE(P’), Fcp, Fcu, GcL, and Ger, need to be solved. However, the points at the 
two sides of the interface between two sub-CEs belong to different SEs. In that case, 
the fluxes across these interfaces may be discontinuous. For this reason, they can be 
calculated by an upwind procedure. For example, 


Fcp = F (uc, ue), (4.40) 


where F refers to any efficient Riemann solver, uçp and ue p denote the values at 
the centroid of the interface FPP’F’, and they are respectively approximated by 


A At 
i Uae omei Poa (4.41) 
A At 
uE = Urne TS tig (4.42) 


As an analogue to the 1D upwind CESE scheme, reconstructed slopes are 
employed in the calculation of the upwind fluxes to eliminate the spurious 
oscillations: 


(ux) rp = ux(A)WBAP™(1, Op» olp) (4.43) 


x 


(ut) pp = ¥x(B)WBAP™(1, Okp, OR) (4.44) 


x 


where the parameters in WBAP limiter are 


pis = (Urp — ULp)/(Ax/2) 92. = ux(B) 
aa ux(A) > Ay 
— Urn — Urp)/(Ax/2) > ux (A) 


ol, = = . 
=, ux(B) ai) 


(4.45) 
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The temporal derivatives are hereafter calculated by the chain rule and Eq. (4.1), 


(ur) pp =-fy (Urp, (už) zp) = 8&y(Uxp, uy(A)) (4.46) 


(U7) ro = ~f: (Uro. (už) rp) — 8y(Uro, uy(B)) (4.47) 


4.2 CESE Schemes on Unstructured Meshes 


In order to extend the capability of the CESE scheme to solve problems with complex 
spatial geometries, CESE schemes on unstructured meshes have been proposed and 
developed. Notably, both central and upwind versions have been developed for hybrid 
meshes of triangular and quadrilateral elements. Hereafter, the second order CESE 
scheme on hybrid meshes [7] is presented here. 


4.2.1 a-a CESE Scheme 


Consider a 2D mesh consisting of quadrilateral and triangular elements as sketched 
in Fig. 4.4. The V; denotes the vertices of the cells. Initially, the value of u, and 
its spatial derivatives, uy and uy, are provided at these vertices. For each cell, C,, 
represents its centroid, and O,, is the mid-point of the cell edge. Meanwhile, V;, V,, 
and A denote points at time step tn—-1/2, tn, tn+1/2, respectively. Other space-time 
points are defined similarly. 

The staggered marching strategy is again adopted here. The scheme marches from 
vertices V; to cell centroids C i fromt = t, to t = ty412. The corresponding CE for C”, 
is either a triangular prism or a hexahedron (Fig. 4.4c and d). SE(C. nD is defined as 
four (e.g., SE (C 1) or five (e.g., SE(C. 5)) planes intersecting at the centers. Similarly, 
in the other half marching step, from t = tn—17⁄2 tot = tn, Cm and Onm form a polygonal 
element of which the centroid is denoted as G;, and the scheme marches from C,, 
to the G;. Note that G; does not necessarily coincide with V;. Thus, interpolation 
is required at each marching step. SE(V, ) is defined as M vertical planes and one 
horizontal plane intersecting at V,, where M denotes the number of C,, linked to 
CE(V;). Each CE is composed of a group of hexahedral sub-CEs. As shown in 
Fig. 4.5, On—1Cm Om V;O,,_,C,,O,,V, is the mth sub-CE belonging to CE(V, ), Gin 
and G, are defined as the centroids of the bottom surface Om—1Cm OmV; and top 
surface OCh O,, V, , respectively. 


E i 
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(c) CE(C1) (d) CEC) 


Fig. 4.4 Definition of conservation elements in CESE schemes based on a hybrid mesh 


It should be noted that the marching scheme from t = t, to t = ty41/2 is a particular 
case of that from f = f,_1. tot = t. From t = t, to t = tayın, the projection of a 
CE on the spatial domain is either a quadrilateral or a triangle, while the projection 
in the other step is a polygon with an arbitrary shape. Without loss of generality, 
we present the marching scheme from t = t,_12 to t = t,. By applying the integral 
conservation law Eq. (4.1) to the mth sub-CE, the flux-balancing relation is obtained 


as 
f u do = f u do — f V-ndo, (4.48) 


c0, V] Om—1Cm Om Vi On—1CnCi, Of, Cn Om 04C), 


mm 


+Om-1V; V/O}, ;+OmVi Vi OF, 


fou 


m—-l~m~m "i 
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Fig. 4.5 The mth sub-CE of CE(V;) 


where V = (f, g) is the vector of the fluxes, and n is the unit outward normal vector 
of the surface. In second-order schemes, the integration over a surface is equivalent 
to the average value at the centroid multiplied by the surface area. Therefore, the 
above equation can be rewritten as 


2 


t D) pd Dall t = A 
Un = _ as >X (A yO FO T Ax® GO) F 25n, (ria = Baka): 
l=1 


(4.49) 


where Sm denotes the area of the top surface, and At/2 = ta — ty-1/2. U.S Um; 
(FY, GW), and (F2, Gm’) represent the average values over On 1m Om V, 

Om—1CmOmVi, Om—1Cm (O 0. and Cy,Om O, es respectively. They can be 
approximated by first-order Taylor expansion in the corresponding SEs. F m-1 and 
F m are the fluxes in the normal direction across the interfaces Om—1 V; V; O, 1 and 
Onm Vi V; oe with the neighbouring sub-CEs. Lm-1 and Ln are the lengths of the line 
segments O,,-;V; and O,,V;, respectively. The summation of the flux-balancing 


equations of every sub-CEs is expressed as 


f u do = 3 f (4.50) 


m= lo 
C1 03C3-.-Cu m-1Cm Om Vi, 


mm 


m—1? 


where Oo = Oy. By substituting Eqs. (4.48) and (4.49) into (4.50), the interface 
fluxes between sub-CEs are cancelled, the conservation law on the entire CE(V; ) 
yields 
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2 


M M M 
u(Gi) = > Um Sm — 3 Sy (‘PFL -— axat) XO Sn (450 
m=1 


m=1 m=1 l=1 


Equation (4.51) is the explicit marching scheme for u (G,). In addition, the treat- 
ment of spatial derivatives is similar to that in cartesian meshes. By using the 
continuous assumption at C',, the following relation can be established 


At 
u(Cin) + ui (Cm) = u(C,,) = u(G;) + ux (G})6Xm + uy (G; )ôYm (4.52) 
for each m, where 6x, = x(Cm) — x(G;), and dy = y(Cm) — y(G;). Apply the same 


relation at C 7 4 > together with Eq. (4.52), the two spatial derivatives are ready to be 
solved using Cramer’s rule, 


ÔUm Ym ÔXm  ÔUm 
5Um ô m+ , ÔXm bUm 

u,(G)) = ie) ee u,(G)) = Le ee, (4.53) 

ÔXm  48Ym ÔXm 45m 

ÔXm+1 ÔYm+1 ÔXm+1 ÔYm+1 

where 
At ; 

Uum = U(Cm) + w (Cm) — u(G;). (4.54) 


Totally, M pairs of derivatives can be calculated in this approach. Shen et al. [7] 
suggested using the weighted average function to obtain the final derivatives as 


M M 
D wW™y™(G') D we” (G') 
u,(G) = = —, u6) = = À (4.55) 
Y we) Y wm 
m=1 m=1 
where 
< (i) 2 (i) a) 
w=! JI [us (G);| + [us (G);| (4.56) 


i=l, im 


Here, the adjustable variable æ (a = 0 ~ 2) controls the dissipation near the discon- 
tinuities. Smaller œ results in lower dissipation, and the function will reduce to an 
arithmetic average when a = 0. Finally, the value at V; is interpolated from G; as 


u(V;) = u(G;) + ux(G))[x(V/) — x(G] + uy (G)[y(V) — yGH]. 4-57 
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The a—a CESE scheme on the hybrid mesh formulated by Eqs. (4.51) and (4.57), 
while robust, is nevertheless sensitive to the Courant number. The following section 
addresses the further improvements to make the scheme Courant number insensitive, 
including the CNI scheme and upwind scheme. These two updated schemes, being 
closely analogous to those presented in Sects. 4.1.2 and 4.1.3, possess the identical 
time-marching scheme for u as provided in the a-a CESE scheme on a hybrid mesh. 
Only the marching scheme for spatial derivatives is modified accordingly. 


4.2.2 CNI CESE Scheme 


As a straightforward extension of the 1D CNI scheme, Shen and Parsani [8] extended 
the CNI scheme to a hybrid mesh. Recall that the central idea of the CNI scheme is 
that when the local Courant number (v) approaches the global Courant number (vo), 
the scheme becomes the a-æ scheme. When v approaches 0, the scheme reduces to 
the non-dissipative counterpart. As depicted in Fig. 4.6, a new point I, is defined 
such that 


x(I) = TAC + (1 = ~) Gn) (4.58) 


yl) = —y(Ci,) + (1 = ZJ). (4.59) 
Vo Vo 


The value of J,, can be obtained using the expansion from G; as 


u(n) = u(x a) = x(G), Yn) = IG), O) gy (4.60) 
Fig. 4.6 Definition of points C x 
in instructed mesh for CNI t ; 
scheme 0; oo" a. O, 


` 
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where u(J,,) can be calculated as 


A 
u(n) = u(x = x (Cm), Yn) — (Cm), >) (4.61) 


Cn 


With the aid of Eqs. (4.60) and (4.61) to each J,,. Every two neighboring T, 
explicitly provide one set of ux (G;) and u (CG). Then, the optimal derivatives can 


be weightedly averaged from these M sets of derivatives. 


4.2.3 Upwind CESE Scheme 


In the upwind CESE scheme, the definition of SE is modified and SE(V;) is the 
polyhedron C\O,C)...Cy OyC; PO; C3...Cy Oy (Fig. 4.4b). The values of u, f, 
and g are continuous within each SE. However, the interface flux Fn between two 
neighboring sub-CEs may be discontinuous. The flux Pm can be evaluated in an 
upwind manner in the normal direction as 


Ên = Ê (uz, ur) (4.62) 


where Ê denotes the Riemann solver, and the values at the center of the interface, 
uz and up, are reconstructed from SE(Cm) and SE(Cm+1) with a similar procedure as 
the 2D upwind CESE scheme on a uniform mesh. The WBAP limiter (Eq. (3.42)) 
is used to reconstruct the derivatives. Once the fluxes across inner interfaces are 
determined, the averaged value on the top surface of mth sub-CE U, in Eq. (4.49) 
can be explicitly computed. We can therefore relate U,, to u(G,) using the Taylor 
expansion as 


Un = u(G;) + ux(G))[x(G,,) — x(G)] + uy (G)[WG,,) — 1G] (4.63) 
Two linear equations can be derived for each set of m and m + 1. Accordingly, a 


similar procedure for averaging the derivatives in the central scheme can be directly 
applied here. 
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4.3 Numerical Examples 


Two canonical problems are selected to demonstrate the performance of the a—a, 
CNI, and upwind CESE schemes on either structured or unstructured meshes. For 
Euler equations, the above schemes can be directly applied to solve each component 
of the conserved variables. The interface fluxes in the upwind CESE scheme can be 
solved as a Riemann problem, where Riemann solvers including Harten, Lax and 
van Leer (HLL) Riemann solver, the contact discontinuity restoring HLLC Riemann 
solver, and the Roe Riemann solver, can be applied to solve the local Riemann 
problem. Since the upwind direction may not be perpendicular to the grid-aligned 
direction, here we use the rotated HLLC Riemann solver. The first example is the 
Mach 3 wind tunnel problem with step. The computational domain is [0, 3] x [0, 
1] with a step of height 0.2 placed at 0.6 from the entrance. The inflow boundary 
condition is imposed at the left boundary, and supersonic outflow is applied at the 
right boundary. The upper and lower boundaries are treated as reflective boundaries. 
Initially, the primitive variables (0, u, v, p) are set as (1.4, 3, 0, 1). Figure 4.7 
depicted the simulations using different CESE schemes with 1800 x 600 structured 
cartesian meshes. Both the CNI scheme and the upwind scheme can capture the shear 
layer instability that is missed in the result from a—a scheme. 

The second example is the double Mach reflection problem. A Mach 10 shock 
propagates to the right and reflects over a solid wedge with an incline angle of 30 
°. A supersonic inflow boundary condition was implemented on the left boundary. 
Reflective and non-reflective boundary conditions were applied on the lower and 
upper boundary, respectively. The problem was solved using unstructured quadrilat- 
eral meshes with a mesh size of 1/400. Figure 4.8 compares the results by applying 
different schemes at computation time t = 0.2. All the three schemes can accurately 
describe the structures while the upwind scheme surpasses with finer resolved feature 
near the wall-adjacent jet. 


bZ 


(a) a-aCESE 


(b) Rotated HLLC CESE 


(c) CNI CESE 
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Fig. 4.8 Double Mach 
reflection problem computed 
by different CESE schemes 
with unstructured 
quadrilateral meshes: density 
contours at t = 0.2 
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Chapter 5 A) 
High-Order CESE Schemes ciecie; 


This chapter is dedicated to the description of high-order CESE schemes. In the 
second-order CESE schemes, the first-order Taylor expansion was employed to 
approximate the unknowns and fluxes within the solution elements. The accuracy of 
the scheme mainly depends on the approximations on the surfaces of the conser- 
vation elements. Analogously, the high-order CESE schemes with Mth-order in 
space and time are generally derived from (M—1)th-order Taylor expansions in the 
solution elements. The high-order CESE schemes use a highly compact stencil. 
Furthermore, spatial and temporal high-order accuracy can be achieved simultane- 
ously. We shall start with constructing a 1D high-order scheme and then extend it to 
multi-dimensional schemes. 


5.1 Construction of High-Order CESE Schemes 


Several attempts have been made to obtain higher-order accuracy for CESE schemes. 
One of the primary advantages of the high-order CESE schemes is the usage of the 
most compact stencil. High-order accuracy is achieved by the approximation with 
high-order Taylor expansions. For example, a(3) scheme [1] has 4th order of accuracy 
by applying second-order Taylor expansion, and it is stable for v < 0.5. Furthermore, 
a(4) scheme [2] with 4th to 5th order of accuracy was developed by defining more 
CEs at each grid point. The CFL number needs to be constrained below 1/3. The 
constructions of high-order schemes from the above two approaches were limited to 
only 1D scenarios. Alternatively, Chang [3] proposed a high-order scheme that can 
be extended to arbitrary order with CFL limited below 1. In this scheme, the even 
derivatives are advanced in the same manner as the conserved variables, whereas 
the odd derivatives are advanced through finite difference. This method has been 
extended to solve Euler equations [4]. In this section, we follow the approaches of Liu 
and Wang [5] in deriving the high-order 1D scheme, Wang et al. [6] for 2D high-order 
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schemes on uniform meshes, and Shen et al. [7] for the 2D high-order scheme for 
hybrid meshes. In these schemes, only the conserved variables are computed through 
the integration on the CE surfaces, while all the physical derivatives are computed 
through the finite difference of lower-order derivatives. In the following text, we 
will present the constructions of the third-order CESE schemes. The construction of 
arbitrary order schemes could be achieved by implementing the same strategy, and 
details can be accessed in Shen et al. [7] and Yang et al. [8]. 


5.1.1 Construction of a Third-Order 1D CESE Scheme 


Recall the 1D scalar problem 


ðu af(u) _ 
ao ee (5.1) 


with f = au. To formulate a third-order scheme, u(x, t) and f (x, t) in any (x, t) € SEG, 
n) are approximated by the second-order Taylor expansion as 


U(x, E) = U} + (Ui — xj) + (UAC = th) + (Ware = xy (t= th) 


1 1 
HUA + suit), Ot) € SEG 62) 


fa DSF +U e-t taa 
1 1 
HSSE + DGE h) œ, t) € (SE); (5.3) 
where tyt, Uxx, Uns f xt» f xx, and f y are the second-order derivatives of u and f, respec- 


tively. These derivatives are assumed constant within the SE. With the aid of the Eqs. 
(5.1) and (5.2), then one has 


(u)! = aus), (Ux)! = alus)", (Un)! = alun)" 64 


It implies that u, uy, Uxx are the only independent variables associated with each 
grid point. Substituting Eqs. (5.2) and (5.3) into the integral form of Eq. (5.1). The 
following algebraic relation can be derived as 


3 


Wj AX + xx) y = Uzr + Ur + Fr — Fr, (5.5) 


where 
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—1/2 2 —1/2 3 —1/2 
Uy a a Ae eye 
Up = & n—1/2 Ax? n—1/2 Ax? n—1/2 
R= THa -3 Udizyp + “ag (Hax) jsp (5.6) 
At pl At n— Ar n j 
F= Ffa tE i + ae i 
At pn- A n A ie 
FR = F j+1/2 + FZ D2 + ag (Sit) ji 


In Eqs. (5.5) and (5.6), the values at t = t,_1/2 are already known. It is necessary to 
calculate (uxx) att =f, before explicitly computing u’. The second-order derivative, 
(Uxx)? , is approximated by a central difference as 


ux) = (Ux)i_1/2 


Ax (5.7) 


(uxx) = 
Then u’; can be directly computed from Eq. (5.5), and (ux); can be calculated 


in the same manner as the second-order a-a schemes using the weighted average 
function. 


5.1.2 Construction of a Third-Order 2D CESE Scheme 
on Uniform Mesh 


Recall the 2D scalar hyperbolic conservation law 


ðu of og 
—+—+—=0. 5.8 
ot ox ody ee) 


which can be cast into the integral form 


$ h-ndS=0 (5.9) 
S(CE(P’)) 


where h = (f, g, u) is the space-time flux vector, n is the unit outward normal vector 
on the surface of the control volume. Inspired by the construction of the high-order 
1D scheme, second-order Taylor expansion is utilized inside the SE to estimate u, f, 
and g. The definitions of SE and CE are kept the same as the second-order schemes 
(Fig. 4.1). For a point (x, y, t) that belongs to SE(P’), 


u(x, y,t) =u(dx, dy, ôt)p, (x, y,t) € SE(P’) (5.10) 
f(x, y, t) = f (6x, dy, 5t)p, (x, y, t) € SE(P’) (5.11) 


g(x, y,t) = g(6x, dy, dt)p, (x,y, tHE SE(P’) (5.12) 
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where X (ôx, dy, ôt)y denotes the second-order Taylor expansion about point N as 
X (x, dy, St)y = Xy + (Xx)ybx + (Xy) Sy + (Xp) dt 


1 2 1 2 1 2 
+ z Xr) x) + 5 (X) O) F z Xn êt) 
+ (Xxy) ySxdy + (Xyr) ySVSt + (Xxr)yôxôt. (5.13) 


and 
ôx =x — xy, ôy = y — Yyy, ôt =t — ty. (5.14) 


Substituting Eqs. (5.10)-(5.12) into (5.8), one obtains 


ur = — fx — By 

Ui = = far = By (5.15) 
Uxt = — fxx — &xy 

Uyt = — fry — Byy 


Together with the chain rule, to compute the derivatives in Eq. (5.13), the only 
unknowns are u, Ux, Uy, Uxx, Uyy, and uyy. By integrating over the surface of the 
CE(P’), Eq. (5.9) leads to 


f n-nas= f h-nds+ f hens + f h-ndS 
S(V) A'B'C'D! ABCD ABBA! 


+f h-nas+ f h-nas+ f h-ndS = 0. 
BCC'B’ CDD'C’ DAA'D’ 
(5.16) 


Thereafter, with the aid of Eqs. (5.10) -(5.12) and (5.16) can be rearranged as 


(Ax)? (Ay)? 1/_ At. At 
U) py +~ Uxx) p + Sa (My) » = (7+ TT Sa) (5.17) 


A Ax A 
-20) +i(= -*.0) (5.18) 
C D 


Ay At -(. Ay At 
0, ’ f 0, , 
4° 4), 4° 4); 
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Ay At A Ay At 
0, + f{0,-—,—] , (5.19) 
c D 


4 

P Ax At ~{ Ax At 

a( , 0, ) e( , 0, ) : (5.20) 
Cc D 


and 
a(dx, 8y, dt) y = (Wy + (Ux) dx + (Uy) yy + (Uy) St 


1 1 
+ z Crx) (8x)? + (Uyy) y OV? + (try) yxy (5.2D 


f (8x, dy, 8t)y = (fn + (F nôx + (fy) yd + ft 


1 1 
+ AUGA + zN Bt)” + (fyt) ySySt, (5.22) 


8 (5x, Sy, St) y = (8)y + (Bx) ydx + (gy) SY + (gi) nôt 


1 1 
+ gee Ox)” + g Enn ON? + (xr) ydx5t. (5.23) 


Prior to explicitly computing (uw), from Eq. (5.17), (xx), and (uyy) y are 
required to be obtained from the finite difference of the interpolations from the 
previous time step. The second-order derivative (uxx), can be estimated as 


(ux)g + a (Uxt)g — Uxda a (Uxt) A 


(uxx) y = = (5.24) 
Uy g At Uy Ux At Uy 

Cae Uxde + FC pe po -~ 7. Od (5.25) 

(Uxx),/ is then computed from simple average as 
Coe ee + (uxx) y 
(Ux) y = —— = (5.26) 
Similarly, the derivative (u vy) is computed as 
+ (uy) p + S (uyi)p (uy) 4 S (yr) 4 
(uy), = F , (5.27) 
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a (uy)c F a (uyi)c (uy) g F (uy) z 


(us), = X l (5.28) 
Uy / + u y / 
(uyy) y = P (5.29) 


Equations (5.17), (5.26), and (5.29) form the procedure to get the values of (u) 


(uxx) p’ and (uy) y at the new time level. To complete the time marching, we need 
to further compute ux, uy, and uyy. The weighted average of the sided-estimations 
computes the first-order derivatives as 


wne al Oe) dell 00, e 
Ux), = 7— |u rM a u rY: Ta ~ Ux) py 
P Ax 2) % 2 Je 4 
= 1 At At 
(Ux) y = Ax 2(Ux) py =u 0, 0, 2 i —u 0, 0, DJ 5 


tae 0,0, £ ieee 2 5.32 
wahe (00%), amy) o 


, 
p? 


; (5.30) 


; (5.31) 


2 1 At At 
(uy); =S 2(uy) — u| 0,0, >), — u(o 0, >), , (5.33) 
(Ux) p = W((ux) p (u$) p» æ), (5.34) 
(uy), = wò) (u$) p> a). (5.35) 


Meanwhile, the mixed derivative uxy is computed from 


+ Ux)p+ Ar (uxt)p — Uxda a Ux) 


(usy), = X , (5.36) 
Coe wet F We on FUu)e (6.37) 
Te Us +F Co F Unda (5.38) 
one wet F eimi Fudo (5.39) 


(usy), = z[ (us) + (ux); + TAR + (Hx); (5.40) 
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5.1.3 Construction of a Third-Order 2D CESE Scheme 
on Unstructured Mesh 


Similar to constructing the third-order scheme on a a uniform mesh, second-order 
Taylor expansion is utilized here, and the variables stored at each grid point are still 
U, Ux, Uy, Uxx, Uyy, and uyy. Now we recollect the definitions of SE and CE in the 2D 
CESE scheme on the unstructured mesh as depicted in Fig. 4.4, and the flux balancing 
relation was derived in Eq. (4.48). From the summation of the conservation law for 
all the sub-CEs, we arrive at the flux balancing relation for CE( V,): 


M M 


f udo=9 f udo- i V - ndo. 


CiO, -Ch Oly MELC p On V; Omi BN p Om 0 Ch+Om-1CnCh O! 


m“ m-l 


(5.41) 


The LHS of Eq. (5.41) represents the integration of u on the top surface of CE( V;), 
i.e., 


M 1 u(G;) 
I u do =u(G})- S+ TERMI = u(G})-S+ $. I I a2 6? 
X 


lollo! m=l~ of Glo! 
ci 01 Cy Om Cm OmG; On-1 


1 a?u(G;) A 


vu(G}) 
ee 2a i PEEN 
+> a2 (Sy)? 4 Te (8x) (Sy)dxdy, (5.42) 


where ôx = x — x(G,), ôy =y— y(G;), and S represents the area of the polygon 
C, 0, ...Cy O'y. The second term in RHS (TERM!) of Eq. (5.42) denotes the inte- 
gration of the high-order terms of u. The two terms in RHS of Eq. (5.41) can be 
respectively expressed as 


M M 
a a 
> [| u do = TERM2 = J` I u(Cm) 4 KEM sw ; HOM (sy) 
x y 

™=1Cm OmV; Om—1 ™=1Cm Om Vj Om] 

1 07u(Cm) 9, 1 82u(Cm) 2 , B2u(Cn) 

pe $x)? + = dy)? 4 8x)(Sy)dxdy, 4 

ae Vti p O + ag OOO (5.43) 


where ôx = x — x (Cm), dy = y — y (Cm). 


“M M 
> f V . ndo = TERM3 = > [FLUX + FLUXP |, 


m=1 m=1 


Cm Om Om Cy + Om—1 Cm Cnn O, 


mm m“ m-l 


(5.44) 


with 


64 5 High-Order CESE Schemes 


2-a2—-a-b ga+b+c c 
m) sm) 1 a Sf (Cm) (m) (m) b| At 
FLUX Say Ax Ay = 
1 > 2 E (a+b Dalbe +D!  axtaybare [ax | [ Xi ] 2 


a=0b=0 c=0 
2—a 2—a—b a+b+ c 
Am) 1 a “g(Cm) (m) an) JÉ At 
Sax Ax Ay awl Be 
Eps 2 (a@+b+Dalbl(e+ D! ax4aybare “1 Fisi | 2 
a=0b=0 c=0 
(5.45) 
2—a 2—a—b +b+e c 
(m) _ At (on) 1 ae f(Cm) (m)]@ an) J| At 
FLUX > AY: Ax Ay. 
2 aa 2 G +b+ Dalbe + 1)! ax@ayPare | “2 I'l ”2 ] 2 
a=0b=0 c=0 
2—a 2—a—b a+b+c c 
fm) 1 a g(Cm) (m)]@ (m)]?[ At 
Sax A Ay. , 
2 D GiPEDaREED! Dar Ha] par] 5 
a=0b=0 c=0 E 
(5.46) 


in which Ax” = = x(Om) — x (Cm), Ay” = = y(On) — y(Cm), Ax = = x(Cm) — 
x(Om-1), Ay! = = y(Cn) — y(Om—1). FLUX | and FLUX", represent the fluxes 
across the adjacent CE surfaces on which the points are defined within SE(C,,). 

By substituting Eqs. (5.42)-(5.44) into (5.41) and moving the high-order terms 
to RHS. Consequently, the value of u(G;) can be expressed by 


1 
u(Gi) = y (TERM2 — TERM3 — TERMI). (5.47) 


Two remarks should be noted regarding Eq. (5.47). First of all, since u, f, and g 
are nonlinearly distributed on the surfaces of CE, the integral terms appeared in Eqs. 
(5.42) and (5.43) are not simple multiplications of the value at the centroid and the 
surface area as what has been implemented in the second-order scheme. In contrast, 
they are now calculated using the coordinate transformation method: 


4 
x=) uN, n), (5.48) 


j=l 


4 
y=) VN: n), (5.49) 


i=1 


where 


1 1 
Nya] aha); Nz = -0 +E) = n), 


4 
1 1 
N3 = g0 TAn), N, = g0 78A tn). (5.50) 


Thereafter, the integration of an arbitrary function $(x, y) over the irregular 
quadrangle A is transformed into the integration over a rectangle (Fig. 5.1) as 
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(X45 Y4) (-1, 1) (1, 1) 


(x3, y3) 


(x yı) (-1, -1) (1, -1) 
(a) (b) 


Fig. 5.1 Sketch of the coordinate transformation 


1 1 
f| 0 as = f an f owe. n), yE, n)]JdE, (5.51) 


-1 -1 


where J is the Jacobian matrix of the coordinate transformation, 


_ ôx, y) a 1 Xi JE 7 Dia Yi ve (5.52) 


4 ON; 
aE, n) Ji 1 my ye 1 Yi On an 


In addition, TERM2 and TERM3 in Eq. (5.47) can be computed from the infor- 
mation at tf = f,_12. TERM! consists of high-order terms (uxx, Uyy, and uxy) of u at 
t = tn. As a result, the high-order derivatives must be solved prior to computing u. 
The idea to calculate the high-order derivatives are analogue to that in Sects. 5.1.1 
and 5.1.2. For example, from the Taylor expansion C i we have 


ux (Ci) — ux (Cn) = uxx (Cy mC -1)-*(Ch)] + Uxy( Ch) D-yCh, )|. (5.53) 


tx (Cine) = Woe (C) = tex (Cin) [Chg = Cin) ] + ery (Ch [EnD — Cin) (5.54) 


On the other hand, Ux(C,,_4)> ux(C,,)s ux (Chg) can be interpolated from the 
previous time level, 


, At 
ux (C4) = Ux (Cm-1) + Uxx (Cm 
i At 
ux (Ch) = Ux (Cn) + Uxx (Cm > 
At 
Ux (Chy) = Us (Cm41) + Uxx (Cm1) -5 (5.55) 


By substituting Eq. (5. 55) into Eqs. (5.53) and (5.54), a pair of equations with two 


unknowns (uxx (C, ), Uxy(C D) are ready to be solved by Cramer’s rule. A simple 
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average (or weighted average to suppress oscillations near discontinuities) can be 
applied for M pairs of derivatives. Once all the second-order derivatives are obtained 
under the same routine, u (G,) can be immediately calculated from Eq. (5.47). The last 
unknowns remain to be computed are uy and uy, which can be calculated similarly 
in Sect. 4.2, e.g., 


u(G‘) + ux(Gi)5x + uy(Gi)dy + Sux. (Gi) (8x)? + duyy (G) 0)? (5.56) 

+uxy (Gi) (dx) (8y) = u(Cm) + ur (Cm) ¥ + Fua (Cm) (2) , l 
where ôx = x(Cm) —x(G;), dy = y(Cm)— y(G,). Similarly, on applying Eq. (5.56) 
to C41, a pair of u and u, can be easily determined. Then, we apply the weighted 
average function to obtain the optimal derivatives. Finally, u(V,) and its derivatives 
are interpolated by Taylor expansion from the information at u(G,). 


5.2 Numerical Examples 


The Shu-Osher problem consists of a Mach 3 shock interacting with an entropy 
wave, which requires the capability of capturing fine structures in the flow. The 
computational domain [0, 10] is discretized by a uniform grid size of 0.05, the non- 
reflection boundary condition is imposed on both sides, and the initial condition is 
described as 


(3.857143, 2.629369, 10.33333) x < 1 


(1+ 0.2sin(5x), 0, 1) x>1 997) 


(p, u, n=] 


The density distributions at t = 7.8 are depicted in (Fig. 5.2). Compared to the 
benchmark solution, which is calculated with a-æ scheme using very high grid reso- 
lution, higher-order CESE schemes are less dissipative and can reveal better flow 
structures. 

The second example is the Mach 3 wind tunnel problem with step. This problem 
has been described in Sect. 4.3. Here, we use the identical numerical setups except 
with a mesh size of 1/250. As depicted in Fig. 5.3, vortex-like structures generated in 
the results using high-order schemes, however, disappeared in the low-order schemes 
due to the increasing dissipation. It is not surprising that the triangular mesh result 
surpasses the quadrilateral ones because the mesh area is only half of the latter when 
the mesh sizes are the same. 
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Fig. 5.2 Density distributions for Shu-Osher problem at t = 1.8. a Entire view, b Enlargement. 
Courtesy of Shen [7] 
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Fig. 5.3 Density contours for Mach 3 wind tunnel with a step at t = 4.0. Courtesy of Shen [7] 
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Chapter 6 R) 
Numerical Features of CESE Schemes E 


In this chapter, some remarks are made on the numerical characteristics of the CESE 
schemes described in foregoing chapters. Due to the special formulation, rigorous 
analysis of the CESE schemes will take more efforts than that of traditional finite 
difference schemes. However, it is possible to extend the widely used modified equa- 
tion analysis, modified wavenumber analysis, and von Neumann stability analysis 
to the CESE schemes. Here, a stability analysis for the second-order upwind CESE 
scheme applied to the linear scalar convection equation is presented. Additionally, a 
rough estimation of the algorithm complexity is presented to show the computational 
efficiency of the CESE method. In order to demonstrate the accuracy of the CESE 
method in comparison with established numerical approaches, simulation results of 
some canonical problems are shown. 


6.1 Stability 


All CESE schemes in this book are explicit time-stepping schemes. Similar to other 
explicit schemes, the numerical stability of CESE schemes is closely related to the 
Courant—Friedrichs—Lewy (CFL) condition. Only when the CFL condition is satis- 
fied, the system can remain numerically stable. Consequently, the time-step size used 
by a CESE scheme has to be restricted according to the upper bound of the stable 
CFL number v for that scheme. 

It has been proven in [1, 2] that the stability conditions for the second-order a 
scheme (Sect. 2.3) and the second-order a—a scheme (Sect. 3.1) are both v < 1. For 
the second-order upwind CESE scheme (Sect. 3.3), the stable upper bound of the 
CFL number is also unity. This stability analysis is first given by Jiang et al. [3] and 
introduced as follows. 

For illustrative purpose, we consider the simplest upwind CESE scheme applied 
to the 1D linear scalar convection equation 
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ə ə 
Z Gog À (6.1) 
ot Ox 


where a is a positive constant. Without using any slope limiter, the upwind CESE 
scheme for Eq. (6.1) can be written in a matrix form as 


n n—1/2 n—1/2 
qj = 0195 i oh On i (6.2) 


with the vector of time-marching variables 


=| ax (6.3) 


and the coefficient matrices (recall Chap. 3) 


If 1i+v 1 — v? I[1—v -1 +? 
=> 3 = ; 6.4 
Q, Al LA pot-Apoe r= z| rarely OA 
where v is the CFL number, which is defined as 
At (6.5) 
v=a—. ; 
Ax 


Following the routine of von Neumann stability analysis, we assume the vector q 
to be a Fourier component with arbitrary wavenumber k such that 


qi =A,(0)e", (=n <60 <x), (6.6) 


where 0 = kAx is the reduced wavenumber, and i is the imaginary unit. Substituting 
Eqs. (6.6) into (6.2) yields 


An (6) = | Q1? + Qpel? A, 12). (6.1) 
Hence, the amplification matrix for each half step is 
M= Oe" 4 Opell? 


O_o @ 2 sð 
Ea i(v* — 1) sins | 


i(l—v)sin$ (2v—1)cos $ +i(v? — 2v) sin $ 


(6.8) 


The spectral radius of matrix M, denoted as p(v, @), is calculated numerically as 
a function of the CFL number (0 < v < 1) and reduced wavenumber (—z < 0 < x). 
In Fig. 6.1, it can be seen that p(v, 0) < 1 holds for all wavenumbers as long as 0 < 
v < 1. Thus, this second-order upwind CESE scheme is numerically stable when v 
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spectral radius of matrix M 


Fig. 6.1 Spectral radius p(v, 0) of matrix M in the von Neumann analysis of the upwind CESE 
scheme. Courtesy of Jiang [3] 


is less than or equal to 1. Also shown in Fig. 6.1, as the CFL number approaches 0 
or 1, the numerical dissipation in the second-order upwind CESE scheme vanishes. 

For high-order CESE schemes, the stability analysis is more complicate and the 
stable upper bound of the CFL number depends on the scheme’s construction. A 
novel approach to construct highly stable high-order CESE schemes was proposed 
by Chang [4], using the same space-time stencil as that in the second-order CESE 
scheme. The same CFL-number constraint for numerical stability (.e., v < 1) is 
retained, which is favourable for explicit time-stepping computations. 


6.2 Efficiency 


The CESE method has a highly compact stencil in space and time, regardless of the 
order of the scheme. This compactness feature provides the CESE method a great 
advantage when parallel computation is performed. 
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For a CESE scheme, the nominal order of accuracy and spatial dimensionality 
will determine the number of independent variables (denoted by K) to be stored 
and updated at each solution point. For a M-th-order accurate (M > 2) and spatially 
D-dimensional (D = 1, 2, or 3) CESE scheme, K can be determined by a precise 
counting and 


M, D 
K(D, M) = M(M +1)/2, D 
M(M + 1)(M + 2)/6 D 


1 
2 (6.9) 
3 


Notably, the number K in the CESE scheme is the same as the number of degrees 
of freedom for each cell in the discontinuous Galerkin (DG) scheme. Therefore, if 
the CESE and DG schemes are of the same order and used to solve the same problem 
on the same mesh, the memory requirements for a CESE scheme and a DG scheme 
are comparable. 

The CESE method requires neither a reconstruction procedure based on a wide 
stencil nor a multi-stage time-integration technique. Besides, no upwind operation 
related to flux calculation is included in central CESE schemes. Nevertheless, the 
CESE method adopts the Cauchy—Kowalewski procedure to express the temporal 
derivatives in terms of spatial derivatives. It also takes time to update the spatial 
derivatives using the algorithms described in Chaps. 3, 4 and 5. The combined effect 
of all these factors determines the computational speed of the CESE method. 


6.3 Accuracy 


Generally, the accuracy of a CESE scheme depends on the order of the Taylor expan- 
sions, i.e., the degree of the approximation polynomials that are used to approx- 
imate unknowns and fluxes within each solution element. To achieve the theoret- 
ical M-th order of accuracy of the scheme in space and time, the (M—1)-th-order 
Taylor expansions are needed. According to practical applications, the second-order 
CESE schemes, including the central and upwind schemes, can usually provide a 
satisfactory trade-off between accuracy and efficiency. 

The ability to capture shock waves and contact discontinuities with high resolu- 
tion is one of the excellent features of the CESE scheme. Here, two canonical flow 
problems are used to demonstrate the accuracy of CESE schemes by comparing the 
numerical simulation results by CESE and other schemes. As the first example, Sod’s 
shock-tube problem [5] is solved by both the CESE method and the finite volume 
method (FVM). The CESE code is an implementation of the second order a—a scheme 
described in Chap. 3, while the FVM counterpart incorporates the second-order 
MUSCL (Monotonic Upstream-centered Scheme for Conservation Laws) recon- 
struction, the van Leer limiter, and the HLLC Riemann solver. A uniform mesh with 
200 cells and a CFL number of 6.9 are employed in both the CESE and the FVM 
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computations. Density and internal energy distributions from x = — 1 tox = l at 
time t = 6.5 computed by CESE and FVM schemes are plotted in Fig. 6.2a and b, 
respectively, in comparison with the exact solution. As shown in Fig. 6.2, a narrower 
contact discontinuity is resolved by the CESE scheme than the FVM counterpart. 
Meanwhile the CESE results match the exact solution better than the FVM results 
at the tail of the rarefaction wave. 


Fig. 6.2 Comparison of 


numerical results of Sod’s X - exact 
2nd-order CESE (200 cells) 
2nd-order FVM (200 cells) 


shock-tube problem by 
CESE and FVM schemes at 
dimensionless time t = 0.5. 
Courtesy of Jiang [3] 


density 


(a) CESE result vs. FVM result: density 


exact 
2nd-order CESE (200 cells) 
- 2nd-order FVM (200 cells) 
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(b) CESE result vs. FVM result: internal energy 
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(b) Upwind CESE result vs. FVM result 


Fig. 6.3 Comparison of FVM and CESE results for Mach 3 step problem: density contours at time 
t = 4.6. Courtesy of Jiang [3] 
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In the second example, the problem of a Mach-3 inviscid flow through a 2D 
tunnel with a forward-facing step inside is solved by three schemes (a—w CESE, 
upwind CESE, and upwind FVM) of second-order accuracy in space and time. The 
initial conditions, boundary conditions, and geometrical parameters are the same as 
that in [6]. The upwind FVM flow solver is provided by Hao et al. [7]. The same 
computational mesh and the same Courant—Friedrichs—Lewy (CFL) number of 6.95 
are employed for all the three schemes. Density fields at time t = 4.0 are shown in 
Fig. 6.3a and b. Compared with the upwind FVM scheme, both the a—w and upwind 
CESE schemes can accurately capture the structure of shock waves, but only the 
upwind CESE scheme can resolve the phenomenon of shear-layer instability (see 
the lower part of Fig. 6.3b), which is not observed in the present a-~ CESE and FVM 
results. 

To assess the accuracy of the high-order CESE schemes described in Chap. 5, 
the Shu—Osher problem and the double Mach-reflection problem were studied inten- 
sively in Ref. [8], where the fourth-order CESE scheme provided a good resolution 
of fine structures in the flows. 
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Chapter 7 R) 
Application: Compressible Multi-fluid ciecie; 
Flows 


Multi-fluid flows involving shock-accelerated inhomogeneities and shock-induced 
instability play essential roles in a wide variety of problems including, but not 
limited to, supersonic combustion [1], inertial confinement fusion [2], and super- 
nova explosion [3]. Numerical simulations of these complex flows prove to be chal- 
lenging in the presence of moving and deformable material interfaces, especially for 
fluids with large differences in their densities or thermodynamic properties. There- 
fore, a discontinuity-capturing, mass-conserving, and positivity-preserving scheme 
is desirable for compressible multi-fluid simulations. 

Qamar et al. [4] implemented the CESE method of Chang [5] and Zhang et al. [6] 
for solving the one- and two-dimensional compressible two-fluid models of Kreeft 
and Koren [7]. Numerical simulations were performed for gas-liquid Riemann prob- 
lems and interactions of air shock waves with inhomogeneities containing lighter 
and heavier gases. In comparison with the non-oscillatory central scheme [5] and 
the kinetic flux-vector splitting scheme [8], the CESE scheme gives better resolution 
of discontinuities. Because of its outstanding ability for discontinuity capturing, the 
studies on Richtmyer—Meshkov instability (RMD) [9, 10] and Rayleigh-Taylor insta- 
bility (RTI) [11, 12] have greatly benefited from the development of this method. 
Those numerical results not only eliminated disturbing factors in experiments (e.g., 
dimensional effect, boundary effect, liquid mist scattering), but also provided more 
ideal working conditions (that cannot be easily realized in experiment) and more 
flow details (that cannot be timely recorded by camera). 

Recently, the upwind time-space CESE method was used extensively in studies of 
shock-accelerated inhomogeneous flows [13-15] and RMI [16-18]. The compress- 
ible two-fluid flows are described by a volume-fraction-based five-equation model 
[19] coupled with the stiffened gas equation of state [20]. Extensive numerical simula- 
tions were carried out using the maximum-principle-satisfying upwind CESE scheme 
[13], which is an improved version of the upwind CESE scheme presented in Chap. 4. 
The maximum-principle satisfying property is achieved by adopting a very simple 
limiter proposed by [21]. Furthermore, the ability to capture contact discontinuities 
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(material interfaces) can be enhanced by employing the HLLC Riemann solver in the 
upwind procedure. As a result, challenging numerical simulations of the gas-gas and 
gas-liquid Riemann problems were successfully performed, in which the mass of 
each fluid component is conserved and the positivity of volume fractions is preserved 
[13]. 

In this chapter, we use extensive numerical examples to indicate that the CESE 
method captures shocks and contact discontinuities sharply without spurious oscilla- 
tions and proves to be a robust and accurate numerical tool for studies of compressible 
multi-fluid flows. 


7.1 Richtmyer—Meshkov Instability 


The Richtmyer—Meshkov instability occurs when an initially perturbed interface 
separating two different fluids is impulsively accelerated. This impulsive acceler- 
ation, in RMI studies especially, is provided by shock waves mostly. Two basic 
mechanisms dominate the RMI development, are the baroclinic vorticity and the 
pressure disturbance. For the baroclinic vorticity, the extent of the misalignment of 
the pressure gradient across the shock with the density gradient across the mate- 
rial interface makes contribution to the perturbation growth. While for the pressure 
disturbance, the wave system plays an important role. 

As the inertial confinement fusion (ICF) cares more about the interaction of a 
converging shock with a disturbed interface, the converging RMI has become an 
imperative [22, 23]. The nature of geometrical convergence in converging RMI, 
however, makes the perturbation development more complicated because of the 
coupling of the Bell—Plesset (BP) effect [24, 25], Rayleigh-Taylor (RT) effect [26, 
27], and the multiple shock impacts therein. 

In perfect agreement with their experimental images [28], Zhai et al. [17] used the 
upwind CESE method to make further discussion about the converging RMI issues. 
Figure 7.1 shows the shock behaviors and interface morphologies after a converging 
shock interacting with perturbed interface. The inner test gas is SF6 and the interface 
initial amplitude is 1 mm. The comparison clearly approved the applicability of the 
method used in simulating converging shock waves and converging RMI issues. 

The Rayleigh-Taylor effect on the phase inversion (RTPI) was the major concern 
in Zhai’s work [17], which is supposed to be helpful in finding a freezing state 
interface in ICF physics. In this work, the influence of the initial amplitude on the 
RTPI was firstly investigated. Figure 7.2 depicts three converging RMI cases with 
different initial amplitudes (Cases I1, I2 and I3 with ap = 1.0, 1.65, and 2.0 mm, 
respectively), according to which, Fig. 7.3 records the time-variation of the interface 
amplitudes. It was found that the amplitude histories can be roughly segmented into 
three stages after the initial shock compression. At the first stage, the amplitude of 
the perturbed interface increases because of the RMI, while during the second stage 
before re-shock, the amplitude reduces owing to the RT stabilization effect which 
is caused by the stronger adverse pressure gradient near the geometry origin. Based 
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Fig. 7.1 The evolution of a single-mode interface with initial amplitude 1 mm. upper: experimental 
schlieren image; lower: numerical schlieren image at the same instants. (IS = incident shock, II = 
initial interface, RS = reflected shock, TS = transmitted shock, SI = shocked interface, RTS = 
reflected transmitted shock, STS = secondary transmitted shock). Courtesy of Z. G. Zhai [17] 


on the amplitude histories, they concluded that whether a normal phase inversion 
occurs or not depends on the competition between the RT stabilization effect and 
RM instability caused by the re-shock. Especially, for a special initial amplitude, 
there is a critical state with a zero-amplitude of the interface at the re-shock. 

Furthermore, a parametric study was performed to evaluate the influences of 
Atwood number, shock Mach number, and initial radius of the interface on the vari- 
ation of critical amplitude with the azimuthal mode number. Their thorough inves- 
tigation well explained the reason why this sensitive RTPI phenomenon was not 
observed in previous converging RMI studies. 

The even more complicated converging RMI at a dual-mode interface was numer- 
ically studied by Zhou et al. [18] using the same upwind CESE methodology. It is 
not surprising that the dual-mode or multi-mode interface is far more complicated 
than the single-mode converging counterpart due to the existence of mode interfer- 
ences such as harmonic generation and bubble merger. By comparing the detailed 
processes of the interface deformation and wave propagation for different dual-mode 
cases together with the development of a referencing single-mode interface, signif- 
icant influence of the phase difference between two basic waves on the instability 
development can be clearly distinguished. While after the re-shock, the discrepancy 
becomes much smaller due to the weak dependence of the re-shocked RMI on the 
pre-re-shock state. 


80 7 Application: Compressible Multi-fluid Flows 


Fig. 7.2 Schlieren images Case |1 Case |2 Case |3 


showing the evolution of the 
single-mode interface. (PI = r, 


SN 


phase inversion, RTPI = ; A 
Rayleigh-Taylor phase AE 
inversion). Courtesy of Z. G. 
Zhai [17] 
864 us 
984 ys 
1356 ps 1357 ps 1355 ps 
1434 ps 1434 ps 1432 ps 
ae 
A 


Critical state 
1468 us 


1623 ps 1623 ps 1623 ps 


The study on dynamics of the dual-mode converging RMI clearly reveals the mode 
coupling effect, in which the growth of the first mode was found to be inhibited, 
promoted, or not influenced, depending upon the first mode amplitude as well as 
the phase difference between the two basic waves. These findings, including all the 
other parametric studies in their work [18], were considered to be of great help for 
designing an optimal structure of the ICF capsule. 
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Fig. 7.3 Time-variation of interface amplitude in three perturbed cases. (PI = phase inversion, 
RTPI = Rayleigh-Taylor phase inversion). Courtesy of Z. G. Zhai [17] 


7.2 Rayleigh-Taylor Instability 


Other than the shock-induced interface instability, the rarefaction-induced insta- 
bility was performed by [29] using the upwind space-time CESE method. Some 
of their numerical schlieren images of the SF¢/air interface instabilities induced by 
the rarefaction waves are presented in Fig. 7.4, in which the rarefaction waves were 
generated by simulating a diaphragm burst in a SF shock tube. The influence of 
rarefaction wave acting time was demonstrated by varying the distance between the 
initial diaphragm and the gas interface (shown in cases SAI, SA3, and SA5); while 
the influence of the rarefaction wave strength was demonstrated by setting different 
pressure ratios between the two sides of the initial diaphragm (shown in SA3, SA7, 
and SA9). 

Generally, after the rarefaction wave sweeps the SF6/air interface from the down- 
side, RTI is triggered and the perturbation on the interface gradually grows. Different 
from the classical RMI where an impulsive action is exerted, and also different from 
the classical RTI where a continuous body force exists, the rarefaction wave imposes 
acceleration to the interface within a limited time. When the rarefaction wave leaves 
the interface, the baroclinic vorticity dominates the later interface evolution and KHI 
occurs on the interface as the further stretching of spikes and bubbles. 

Based on the numerical results provided by the upwind CESE method, theoretical 
models were modified in Liang’s work [29] considering the time-varying acceleration 
and the growth rate transition from RTI to RMI. It was also found that the interface 
perturbation can be more unstable under the rarefaction wave condition than under 
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Fig. 7.4 Schlieren images of the SF¢/air interface evolution induced by rarefaction waves in cases 
a SA1, b SA3, c SA5, d SA7, and e SA9. Numbers represent the time in ps. f schlieren images 
of the interface evolution in cases SA1 (left), SA3 (middle) and SA5 (right) at dimensionless time 
10.3. (RWH = rarefaction wave head, RWT = rarefaction wave tail). Courtesy of Y. Liang [29] 


the shock wave condition due to the larger amount of vorticity deposited by the 
continuous pressure gradient. 


7.3 Shock Refraction 


The upwind space-time CESE method has also been used in simulating the 
shock refraction phenomenon at an inclined air/helium interface in the cylindrical 
converging shock scenario [16]. Figure 7.5 presents the experimental and numer- 
ical schlieren images of the shock refraction, in which good agreement between 
the two was achieved. Moreover, the numerical results apparently provide much 
cleaner images than their experimental counterparts (especially the region behind 
the deformed gas interface). 

It is observed that, during the incident shock wave converging, when the shock 
velocity increases and the incident angle (with respect to the gas interface) decreases, 
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Fig. 7.5 Sequences of experimental (lower) and numerical (upper) schlieren frames showing the 
evolution of a converging shock wave refracting at a 45° tilted interface. (i = incident shock, m = 
material interface, t = transmitted shock). Courtesy of Z. G. Zhai [16] 


the wave pattern is found to transit from one to another. In the work of Zhai et al. 
[16], two interfaces with 45° and 60° initial incident angles were examined. For the 
45° case, the shock pattern was found to transit from free precursor refraction (FPR) 
to bound precursor refraction (BPR), and then to regular refraction with reflected 
shock (RRR); while for the 60° case, it was found to be from twin von Neumann 
refraction (TNR), to twin regular refraction (TRR), to free precursor von Neumann 
refraction (FNR), and finally to FPR. These clearly depicted transition sequences 
that never occur in the planar shock scenarios, greatly enriched the shock refraction 
classification. 


7.4 Shock-Gas-Bubble Interaction 


Here, we review the numerical studies of shock-accelerated inhomogeneous flows 
conducted by Shen et al. [13], Fan et al. [15], and Guan et al. [14] first. Figure 7.6 
shows a general schematic of the computational setting of a shock-accelerated inho- 
mogeneity. The homogeneity can be gas bubble [13], water column [13] or droplet 
[14], and the shape of the initial homogeneity can be circular, square, rectangular, 
and even triangular [15]. 
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Fig. 7.6 Schematic of the initial setup for the interaction of an incident shock wave and a gas/water 
inhomogeneity 


In the shock—-helium-cylinder case [13], a planar shock propagates from left to 
right at Mach 1.22 and impacts the helium cylinder with a radius of 25 mm. The 
cylinder morphology after shock impact is shown in Fig. 7.7, in which the gas inter- 
face is sharply captured. The vortices, deposited due to baroclinicity, are clearly seen 
developing on the interface. As is shown, a salient feature of the helium cylinder 
morphology is the piercing of a middle air jet that finally separates the whole cylinder 
into two symmetric lobes. 

A quantitative validation of the CESE method in simulating the shock—helium- 
cylinder interaction was performed by comparing the trajectories of the upstream 
point, the jet point, and the downstream point of the helium cylinder with the ones 
obtained by the level set method [30]. Good agreement of the results was achieved 
between the two methods, as depicted in Fig. 7.8. 

The maximum-principle-satisfying upwind CESE scheme was applied by Fan 
et al. [15] to present a comprehensive study of jet-formation phenomenon in the 
interaction of a planar shock with a variety of heavy gas inhomogeneities. Comparing 
with the above light helium cylinder scenario, the heavy R22 cylinder experiences 
an obvious different deformation pattern. Figure 7.9 shows the cylinder morphology 
of the R22 cylinder with a radius of 25 mm impacted by a planar incident shock with 
Mach number 1.22 (to numerically reproduce the experiment conducted by Haas 
and Sturtevant [31]). The shocks that transmitted into this heavy cylinder were found 
converging at the downstream pole. The ensuing shock pattern and pressure gradient 
leads to a jet at the downstream cylinder surface. 

Fan et al. [15] summarized that the heavy gas jet forms not only in the scenario 
of circular cylinder, but also in a series of differently shaped heavy gas cylinders, 
including square, rectangle, and triangle. Fan’s square case reproduced the experi- 
ment performed by Luo et al. [32] where the shock Mach was 1.17 in air and the side 
length of the SFe square was 56.6 mm. To reveal more details of the shock converging 
patterns in a small region, a slightly over-fined mesh was used in their simulation, 
which greatly reduced the numerical viscosity, making the secondary vortices at the 
late stage more pronounced than its corresponding experimental images. Other than 
that, numerical results shown in the lower halves of Fig. 7.10 agree quite well with 
the upper halves’ experimental images, including the gas interface morphology, the 
shock patterns, and the consequent jet formation at the middle leeward surface of 


7.4 Shock—Gas-Bubble Interaction 85 


0.03 |- 0.03 F 
0.02 F 0.02 
0.01} 0.01 

OF 0 
-0.01 -0.01 
-0.02 + 0.02 }- 
-0.03 f A i 0.03 etry 

0 0.01 0.02 0.03 0.04 0.03 0.05 0.06 

(a) time = 199.5 us (b) time = 399.7 us 

0.03 0.03 
0.02 0.02 
0.01 0.01 F 

0 Or 
0.01 -0.01 F- 
-0.02 0.02 - 
0.03F, | 1 ! 00S FS 1 1 1 iit 

0.05 0.06 0.07 0.08 0.09 0.07 0.08 0.09 0.1 0.11 0.12 
(c) time = 600.6 us (d) time = 800.2 us 


1 L ri ee E Lo 
0.1 0.12 0.14 0.14 0.16 


(e) time = 999.5 us (f) time = 1200.0 us 


Fig. 7.7 Deformation history of a helium bubble impacted by a Mach 1.22 shock. Courtesy of H. 
Shen [13] 
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Fig. 7.8 Position history of the upstream point, the jet point, and the downstream point of the 
helium bubble. Courtesy of H. Shen [13] 
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Fig. 7.9 Shadowgraphs of the evolution of an R22 circular cylinder under the impact of an incident 
shock with Mach number 1.22. (is = incident shock; dts = direct transmitted shock; rs = reflected 
shock; rts = reversed transmitted shock; s/ = slip line.) Courtesy of E. Fan [15] 


the square. Based on the detailed information provided by the CESE simulation, 
the mechanism of jet formation was revealed. It was found that for all the cases 
they simulated, the key factor in forming a downstream jet is the formation of type II 
shock-shock interaction [33]. According to this, a geometrical criterion was proposed 
to determine whether a jet will be formed [15]. 


7.5 Shock—Water-Droplet Interaction 


The ability of the upwind CESE method to handle gas-liquid interfaces was initially 
demonstrated in Shen’s work [13]. The setup of this problem resembles that of 
the shock—bubble interaction (see Fig. 7.6). Due to the high-density ratio and large 
difference between the thermodynamic properties of air and water, more challenges 
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Fig. 7.10 Comparison of experimental images (upper) by Luo et al. [32] and the numerical results 
(lower) of the planar incident shock impact on a square cylinder using CESE method. (is = incident 
shock, rs = refracted shock, ts = transmitted shock, ms = Mach stem, gi = gas interface, dts = 
direct transmitted shock, lts = lateral transmitted shock, ds = diffracted shock, uvr = upstream 
vortex pair, dvr = downstream vortex pair). Courtesy of E. Fan [15] 


are there in simulating the shock—water-column interaction. In this case, the radius of 
the water column was 2.4 mm anda shock wave with Mach number 1.47 was launched 
to sweep over the column from left to right. Upon the interaction of a shock with water 
column (or water droplet), two essential issues are extensively discussed. One issue 
concerns the immediate interaction of the shock wave (or rarefaction wave) with the 
water column (or droplet). It covers an extremely short period of time compared to 
the whole water column breakup, in spite of its vital role in the following interface 
deformation. The other issue concerns the long-period water column deformation, 
which covers from the initial protrusion growth on the water column surface to the 
later large interface distortion and mass loss of the bulk of the column. At the early 
stage, as shown in Fig. 7.11, the shock propagates through the water column as 
though passing over a rigid cylinder, nearly no interface deformation is detected. 
However, the shock waves as well as rarefaction waves bounce back and forth inside 
the water column, making the internal pressure change dramatically and becomes 
highly heterogeneous. After the incident shock has passed, as shown in Fig. 7.12, 
atomized water is gradually stripped away by the high-speed post-shocked flow. The 
water column develops into a crescent shape and the downstream was covered up by 
the transversely spreading atomized water. Note that the instabilities were captured 
in detail in this CESE simulation. The history of drag coefficient of the water column 
was also derived from the numerical results, as shown in Fig. 7.13. Again, good 
agreement was obtained in comparison with numerical simulations of Chen [34] and 
experimental data of Igra and Takayama [35]. 
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Fig. 7.11 Sequences of pressure contours (unit: Pa) in the early stage of a shock—water-column 
interaction. Courtesy of H. Shen [13] 


7.5 Shock—Water-Droplet Interaction 89 


oo ove band s 
0.018 coal 001s 
omt A. Ld 3 
0.005 - 0.005 0.005 + J 
0.005 ooo5F j 0.008 Ss 
omp 001 onp 
omst 0015 omst 
0.02 |- oot 0.02 + 
1 4 1 i t 1 4 1 4 L 4 4 
002 rr] oy Oe 002 a an oy, om EA 0.02 Do o M 002 


(a) time = 50.4 us 


002p oxp 
0.015 h osp 
o0 f 001 
~ E 
0.005 f- = = 0005 
> OF > oF 
0.005 - > 2005 
a © 
‘0.01 0.0) 
omst aose 
oot ont ovf- 
zs = 4 ~ L 4 L ds 
o 0.01 g 0.01 00 0 oor oo 0.03 004 0 091 002 0.03 004 
x x x 
(e) time = 400.4 us (Ê) time = 450.0 us 
"j 


a 
ay 


x 7 
7 wy 
o Tor oo Oa To o oor oor 0 Tor o 
(g) time = 500.2 us (h) time = 550.1 us (i) time = 600.5 us 


Fig. 7.12 Evolution of the numerical schlieren plot during the stripping breakup of the water 
column. Courtesy of H. Shen [13] 


The case of a spherical water droplet deformation under planar shock impact was 
studied by Guan et al. [14] using axisymmetric simulation. As shown in Fig. 7.14, the 
initial protrusions emerging on the droplet surface were clearly captured by the CESE 
methodology. The internal flow pattern of the droplet was visualized numerically in 
this study. It is seen in Fig. 7.15 that the internal flow pattern of the droplet is 
established soon after the impact by the incident shock and is held steady for a long 
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Fig. 7.13 Drag coefficient Cg history of a water column impacted by a shock versus dimensionless 
time t“. Courtesy of H. Shen [13] 


S| 


(0, 0) 


Fig. 7.14 Comparison of the numerical (lower part) and experimental (upper part) results at 
different instants. (EQ = equator; WS = windward stagnation; LS = leeward stagnation; C = 
corrugation; KHI = Kelvin-Helmholtz instability, AT = atomization). Courtesy of B. Guan [14] 


time. For the first time, a saddle point was observed in this internal flow pattern. 
Further, a simple theory was proposed to correlate the stationary position of the 
saddle point with the Mach number of the incident shock. 
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Fig. 7.15 Morphologies and internal flow patterns of a water droplet under shock impact with 
different shock Mach numbers Ms. Contours show the density of air with units kg/m°. Courtesy of 
B. Guan [14] 
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Chapter 8 A) 
Application: Detonations ciecie; 


Detonation is a shock-induced combustion in which chemical reactions are closely 
coupled with shock waves. The shock wave compresses the reactant with an abrupt 
increase in temperature and pressure, initiating the reactants to be burnt into products. 
The intense heat release permits the high propagating speed of the shock wave to be 
sustained. It is fundamental research related to both the safety industry and propul- 
sion systems. For most explosive mixtures, detonation wave speeds are formulated 
by Chapman-Jouguet (CJ) theory. Typical detonation velocities for gaseous mixtures 
generally range from 1400 to 3000 m/s. Behind the shock, the time scale for reactions 
is commonly on the order of microseconds or even less. Furthermore, the detona- 
tion front is intrinsically unstable, forming transient multi-dimensional structures. 
Many studies revealed that high resolution is necessary to resolve the essential deto- 
nation structures. Due to its complex nature and multiple time scales, detonation is 
thus a challenging problem for solvers on shock-capturing capability, robustness, 
and computational efficiency. This chapter will present several essential aspects of 
detonation research by applying the CESE schemes. 


8.1 Gaseous Detonations 


In detonation simulations, the requirement for the computational resources is usually 
very high because the satisfactory resolution is necessary to compute the coupling 
between flow and chemical reactions and to resolve the detonation structures. As 
a result, applying detailed chemistry is usually limited to problems with simple 
and small domains or particular situations. Otherwise, simulations with simplified 
models are preferred to gain insight into the physics without loss of credibility. For 
illustrative purposes, we consider a system for ideal gas and a simplified chemical 
model. The 2D inviscid reactive compressible flow equations can be expressed as 
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where the conserved variable vector U, the flux vectors F, G and source term vector 
S are defined as U = [p, pu, pv, E, pd)’, F = [pu, pu? + p, puv, (E + p)u, pau)’, 
G = [pv, puv, pv? + p, (E + pẹ, pav]’, S = [0, 0, 0, 0, òl”. The symbols p, 
u, v, p, E, and à denote density, velocities in x and y directions, pressure, the total 
energy per unit volume, and mass fraction of the reactant, respectively. The perfect 
gas law was adopted here as p = pRT, where R is the gas constant and T is the gas 
temperature. The total energy is 


_. l 2, 2 
E + sp(u? +v?) + pQ, (8.2) 


where y and Q are the specific heat ratio and the specific heat release. The chemical 
reaction rate is formulated by the one-step Arrhenius model as 


ò = —kpheFa/RT (8.3) 


where k is the pre-exponential factor, and E, denotes the activation energy. 


8.1.1 Ignition of Detonations 


Ignition of detonation refers to the formation of a detonation wave through either an 
instantaneous onset or the transition from deflagration. Shen and Parsani [1] studied 
the direct ignition of spherical detonation by using the upwind CESE scheme to 
solve Eq. (8.1). The distance between the shock wave and the position where half 
of the reactant is consumed can be defined as the half-reaction length. A uniform 
resolution with 20 meshes per half-reaction length was adopted, and the number 
of mesh points is 1.6 billion. It is well-known that the dynamics of detonation are 
sensitive to the activation energy. Thus, three activation energies representing stable 
(Ea = 15), mildly unstable (Ea = 27), and highly unstable (Ea = 50) detonations 
were studied. For different E4, k was adjusted to fix the half-reaction length in unit 
length. A hot spot with p, = (y — DE;, Tẹ, = 20.0, and a radius of 1 was used to 
form a blast wave to ignite the mixture. All variables above are nondimensionalized 
with respect to the state of the unburnt reactant. 

The peak pressure history behind the shock for 1D detonations are depicted in 
Fig. 8.1, and typical 2D contours are plotted in Fig. 8.2. For stable case (Fig. 8.1a), 
the subcritical, critical, and supercritical regimes are observed successively when 
E, increases from 1.8 x 10* to 1.0 x 10°. In the critical regime, E, = 2.0 x 104 
(Fig. 8.1b), Psn first decays and then rapidly increases due to the formation and 
amplification of the pressure pulse behind the leading shock. The 2D result shows a 
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Fig. 8.1 Peak pressure history at the shock front (Psn) as a function of position and source energy 
(Es), fora Eg = 15, b Eg = 15, Es = 2 x 10t, c Ea = 27, d Ey = 50. Courtesy of H. Shen [1] 


highly oscillatory pattern for P,, due to the development of multidimensional insta- 
bilities, and the run-up distance deviates from 1D solution. Three critical energies 
were detected for the mildly unstable case (Fig. 8.1c), The authors further observed 
the fourth critical energy when refining the search using a minor incremental Fs. 
However, there might be a unique critical energy in the 2D case. The inconsistency 
between 1 and 2D becomes more apparent for the highly unstable cases. For all the E, 
tested, 1D detonation eventually failed (Fig. 8.1d). The results suggest that, without 
large overdriven, the 1D detonation propagation through auto-ignition is impossible 
for the highly unstable cases. Contrarily, the critical energy is approximately 1.5 x 
10° for 2D cases (Fig. 8.2d). The key factors dominating the direct ignition of 2D 
detonation are the unsteadiness arising from the decay of the leading shock, heat 
release from the chemical reaction, and the inherent multidimensional instabilities. 
In 1D detonations, only the first two aforementioned factors exist. Detonation fails 
when the excessive unsteadiness overtakes the heat release. The transverse waves 
induced by the multidimensional instabilities create local over-driven detonations. 
On the one hand, these local over-driven detonations facilitate the propagation of the 
global detonation; on the other hand, induce stronger expansion waves to increase 
the risk of failure. The competition among these three factors should be considered 
to provide a comprehensive model of direct initiations. 


98 8 Application: Detonations 


2000 2000 


1500 1500 


1000 


1000 


500 500 


0 500 1000 1500 2000 - 
0 500 1000 1500 2000 


(a) (b) 


2000 2000 


1500 


1500 


1000 1000 & 


500 


0 500 1000 1500 2000 0 500 1000 1500 2000 
(c) (d) 


Fig. 8.2 Maximum pressure history in two dimensions. a E4 = 15, E, = 2.5 x 104, b E, = 27, 
E; = 4.5 x 10+, ¢ Ea = 50, Es = 1.0 x 104, d Eg = 50, Ey = 1.5 x 10*. Courtesy of H. Shen [1] 


Another situation of detonation ignition occurs when the gaseous reactant in the 
unconfined space is ignited through the detonation wave emerging from a small 
tube. When the detonation wave diffracts from a tube into the unconfined space, the 
substantial expansion will perturb the detonation front. If tube size is small enough, 
the detonation wave will quench (subcritical), otherwise, the detonation will continue 
to propagate (supercritical). Shi et al. [2] investigated the 2D detonation diffraction 
using the a—æ CESE scheme. Two scenarios were considered. If the detonation wave 
inside the tube does not contain transverse waves, it is classified as the planar case. 
Otherwise, it is classified as the cellular case. In these simulations, Eqs. (8.1)—(8.3) 
were employed with Ea = 24 to assure that 1D detonation wave is stable (For Ea 
higher than 25.27, 1D detonation is intrinsically unstable). A resolution with 24 
meshes per half-reaction length was tested converged for current settings, and the 
total mesh points are around 410 million. A first glance at the results revealed that the 
critical channel width for cellular scenarios is smaller than that for the planar case 
(Table 8.1). Four cases, including supercritical and subcritical for planar and cellular 
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Table 8.1 Critical channel — = 
L limit limit 
width for Eg = 24 ower lim Upper limi 
Planar 100 110 
Cellular 75 85 


scenarios, were carefully discussed by examining wave structures and Lagrangian 
particles that are initially scattered in the flow. 

For planar scenarios, the disturbance originates from the corner propagates 
towards the symmetric line. The diffracted shock wave decelerates and decouples 
with the flame front. In the shocked reactant, a compression wave is formed and 
amplified through the temperature gradient towards the leading shock. The compres- 
sion wave is indicated by the successive abrupt change in temperature for particles 
E7-E9 (Fig. 8.3). When this compression wave merges with the leading shock, deto- 
nation may be facilitated if the strengthened shock wave exceeds a critical value. 
On the other hand, the aforementioned disturbance reflects from the symmetric line, 
and the reflected rarefaction (Fig. 8.4) wave may simultaneously affect the shocked 
reactant and hinder the formation and amplification of the compression wave. For the 
sub-critical case, i.e., w = 100, re-initiation is reproduced if the boundary condition 
is modified such that the rarefaction wave does not reflect towards the leading shock. 
Therefore, for the planar cases, re-initiation is attributed to the competition between 
the coalescence of the amplified compression wave with the leading shock and the 
strength of the reflected rarefaction wave. 

For cellular scenarios, the re-initiation patterns are different from the planar cases. 
Despite the number of transverse waves decreases when the detonation wave expands 
into the unconfined space. It is interesting to note that from Fig. 8.5, a complete 
decoupling is not observed during the propagation of the supercritical case (w = 
85) (a complete decouple is observed before re-initiation in the planar w = 110 
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Fig. 8.3 Temperature traces of particles along the re-initiation direction for planar detonation 
diffraction: a w = 110 and b w = 100. Courtesy of L. S. Shi [2] 
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Fig. 8.4 Numerical schlieren images for planar case with w = 100. Red dash lines indicated the 
head of rarefaction. Courtesy of L. S. Shi [2] 
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Fig. 8.5 Numerical soot foils for a planar detonation diffractions, w = 110, b cellular detonation 
diffractions, w = 85. The red dashed lines indicate the disturbance line predicted by Skews’ model. 
The blue dash-dot line is the re-initiation path along which particles £; are initially located. Courtesy 
of L. S. Shi [2] 


case). Detonation wave sustains through the continuous formation of local overdriven 
detonation when transverse waves collide. Previous experimental research revealed 
that the critical tube diameter for unstable mixtures is de © 13A (A: detonation 
cell size), but de ~ 30A or even more were measured for stable mixtures highly 
diluted with argon. The distinct difference between the re-initiation mechanisms 
planar and cellular detonation diffractions elaborated in current simulation provides 
a plausible explanation of the discrepancy of the correlations on critical tube sizes. In 
unstable detonations, where the formation of local hot spots is essential to sustain the 
detonation, the correlation between the d,. and A is established. It is well-known that 
the cellular patterns in stable detonation are regular, with transverse wave velocity 
approximately being the sound velocity in the burned products. This corresponds to 
the planar scenarios. Detonation re-initiates if the influence of the reflected rarefaction 
wave is negligible. Thus the d, / à correlation may not be necessary. 
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8.1.2 Dynamics of Detonations 


As discussed in Sect. 8.1, multi-dimensional instability is a fundamental phenomenon 
in detonation. In this section, we will elaborate more on the studies on detonation 
propagation in one, two, and three dimensions using CESE schemes. 

Wang et al. [3] have implemented 2D CESE on detonation reflection on a wedge 
using a detailed reaction model. Later, in the study on the propagation modes of 
stoichiometric H2/02 3D cellular detonation in a square duct [4], two types of prop- 
agation modes were verified, which depended on the configuration of the initial 
conditions. Unreacted pockets with complex structures were also captured by using 
the improved CESE scheme. 

Numerous numerical studies based on detailed chemical mechanisms have quanti- 
tatively reproduced the key structures of multi-dimensional detonations. They expe- 
rienced barriers to predict the cell size correctly. For example, for H2/O2 mixture 
highly diluted by argon, the regular cell sizes predicted from simulation using detailed 
chemical mechanisms were approximately half of those measured in experiments. 
One of the important physical phenomena in high-speed flow, the vibrational non- 
equilibrium of molecules, was usually ignored in early numerical studies. Behind the 
detonation wave, the time scales for vibrational relaxation and chemical reaction are 
comparable. Shi et al. [5] re-examined this problem with the consideration of vibra- 
tional relaxation and its coupling with reactions. The convection of vibrational energy 
was integrated into governing equations, and the translational-rotational energy and 
the vibrational energy were separated in the formulation of the total energy, which is 
composed of enthalpy of formation, translational-rotational energy, kinetic energy, 
and vibrational energy. The energy exchange rate was calculated by the Landau-Teller 
model. 

Four scenarios were compared in 1D and 2D simulations. (1) The mixture is 
assumed to be thermodynamically equilibrium. (2) Vibrational relaxation is incor- 
porated, and the translational temperature is kept as the dominant temperature of 
the chemical reactions. (3) Vibrational relaxation is considered, and the geometric 
average temperature is used to compute the reaction rates. (4) The physically consis- 
tent vibration-chemistry-vibration coupling model is adopted to account for the effect 
of vibrational non-equilibrium on chemical reaction rates. Examples of 1D detona- 
tion structures are depicted in Fig. 8.6. Model (2) exhibits a temperature overshot 
behind the shock compared to the equilibrium case because the mixture remains 
vibrational cold. The vibrational temperature approaches the equilibrium state grad- 
ually, and a noticeable disparity still exists when exothermic reactions are triggered. 
Since translational-rotational temperature controls reactions in model (2), the over- 
shot in temperature results in a marginal decrease in the half-reaction length (6). In 
model (3), the geometric averaged temperature is used to calculate the reaction rate. 
The onset of severe chemical reactions occurs only when the vibrational relaxation 
approaches equilibrium. The half-reaction length increases by 1.64 times that in the 
model (1). A similar trend is observed in the case using model (4). The detonation 
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Fig. 8.6 Temperature and Hp distribution with initial condition of 0.1 atm and 300 K. Left: scenarios 
1 and 2. Right: scenarios 1 and 3. Courtesy of L. S. Shi [5] 


cell widths are predicted from 2D simulation using models (1) and (2) are approxi- 
mately the same. The cell widths using models (3) and (4) are approximately 1.32 and 
1.34 times the value of that in the model (1). It reveals that the vibrational-chemical 
coupling effect may be one of the factors responsible for the failure to correctly 
predict the cell size. 

Later, Uy et al. [6] further examine the effect of vibrational non-equilibrium on 
the one-dimensional instabilities using the CESE scheme. In this study, 1D piston- 
supported detonation case with fixed non-dimensional specific heat release Q = 
50, the ratio of specific heats y = 1.2, and characteristic vibrational temperature 
6 = 20. The ratio between the chemical time scale t, and the vibrational time 
scale ty was selected as a control parameter. Park’s two-temperature model was 
used to evaluate the effect of vibrational non-equilibrium on the chemical reac- 
tion rate. The neutral stability limit under the vibrational equilibrium assumption is 
approximately Ea = 26.47, above which the detonation is unstable with amplified 
oscillation. Numerical tests reveal E, = 27 with the thermal equilibrium assump- 
tion being longitudinally unstable. Nevertheless, for finite vibrational relaxation rate 
(i.e., the T” = t,/ty Æ oo), the amplitudes of pressure oscillations decay for smaller 
t“ = 3,5, 7, and the oscillation decays much faster for as t” decreases. The results 
imply that the propagation of detonation is stabilized because of the vibrational non- 
equilibrium. Theoretical analysis using linear stability analysis (LSA) under different 
overdriven factors f =(D/Dc;)’ was also provided in Uy et al. [7]. The neutral stability 
limit predicted using the CESE scheme agrees very well with LSA (Table 8.2). The 
accuracy of CESE in solving detonation problems is further confirmed. 
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Table 8.2 Comparison of the neutral stability limit (NSL) and the period of oscillation (PO) 
computed by LSA and numerical simulation. Case I: f = 1, thermal equilibrium (eq). Case II: 
f = 1, thermal nonequilibrium (neq), t“ = 5. Case III: f = 1, neq, t“ = 7. Case IV: f = 1, neq, t“ 
= 9. Case V: E4 = 50, eq. Case VI: Eg = 50, neq, tT” = 5. Case VII: Eg = 50, neq, tT” = 10. Case 
VII: Ea = 50, neq, t” = 20. Courtesy of C. K. Uy [7] 


Case Linear stability analysis Numerical simulation 
NSL PO NSL PO 

I Eq = 26.46 10.63 Eq = 26.47 10.64 
II Eq = 27.13 12.14 Eq = 27.14 12.15 
il Eq = 26.99 11.77 Eq = 27.02 11.74 
IV Eq = 26.90 11.54 Eq = 26.92 11.51 
Vv f = 1.62 8.02 f = 1.62 8.03 
VI f =1.533 9.51 f =1.554 9.52 
VII f = 1.582 8.80 f = 1.581 8.83 
Vill f = 1.60 8.46 f = 1.598 8.43 


8.1.3 Rotating Detonation Waves 


Because of the rapid compression of the mixture and the quasi-constant volume 
combustion, the high thermodynamic efficiency of detonation makes it favourable 
to be employed in developing potential detonation engines. The rotating detonation 
engine (RDE) is one of the promising candidates in which the detonation wave propa- 
gates circumferentially, and the reactant is continuously injected into the combustion 
chamber. 

Figure 8.7 depicts the flow features of RDE simulated with the upwind CESE 
schemes. The 2D simulation is essentially the unwrapped RDE channel from a 3D 
geometry. The top is the RDE outlet described by a non-reflective boundary. The 
reactant is injected into the combustor from the bottom through micro-Laval nozzles. 
The ideal injection was assumed that all the bottom areas could inject fresh reactant. 
Typical flow features of RDE can be observed, characterized by a rotating detonation 
wave, triangular reactant layer, oblique shock, and slip line. Furthermore, Wang et al. 
[8] used an improved CNI 2D CESE scheme to study the kerosene/air RDE. The 
results suggested that decreasing the inlet/wall ratio would reduce the strength of the 
detonation wave, and the reaction zone would be elongated. With further decreasing 
the air ratio, the burned gas emerges in the triangular zone, and the detonation wave 
fails to be self-sustained. 
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Fig. 8.7 Density contours of 2D and 3D RDE simulations using the upwind CESE scheme 


8.2 Two-Phase Detonations 


Liquid fuels and energetic powders are of advantages of high energy density and easy 
storage, which makes them common in practical propulsion systems. Detonation- 
based propulsion systems using these non-gaseous fuels involve reactive multi-phase 
high-speed flows. Although gaseous, liquid and/or solid phases are involved in these 
detonation phenomena, it is always referred to as two-phase detonation from the view- 
point of reactants, namely gas-liquid two-phase detonation and gas—solid two-phase 
detonation. Compared to gaseous detonations, two-phase detonations are character- 
ized not only by shock waves and fast combustion, but also by multi-phase interac- 
tion and multi-scales, leading to great difficulties in numerical simulations. With the 
successful applications of the CESE method in gaseous detonation simulations, it is 
imperative to extend the CESE method to two-phase detonation problems and test 
its capabilities of accuracy and robustness in two-phase detonation simulations. 

There are mainly two kinds of frameworks that are used to address the discrete 
phase in two-phase detonation simulations, namely the Eulerian—Eulerian framework 
and the Eulerian—Lagrangian framework. Between these two frameworks, the Eule- 
rian—Eulerian framework, which is also referred to as the two-fluid model/method, 
is more common and easier in implementation and extension in two-phase deto- 
nation simulations. The discrete phase (the liquid phase in gas-liquid two-phase 
detonations or the solid phase in gas—solid two-phase detonations) is considered as a 
special continuum, such that continuum mechanics can be employed to describe the 
bulk motion of the discrete phase. Then, the discrete phase becomes another “fluid”, 
and the gaseous-phase flow and the discrete-phase flow can be solved by similar 
approaches. As for the Eulerian—Lagrangian framework, the gaseous phase is solved 
as in the Eulerian—Eulerian framework, but every discrete liquid particle (droplet) or 
solid particle is tracked individually by Newton’s laws of motion. 

In two-phase detonation modelling, it is a general way to assume that the 
particle (liquid or solid) of the discrete phase is small, spherical in shape and 
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uniformly suspended in the gas atmosphere. Additionally, the two-phase suspen- 
sion is assumed to be diluted enough to neglect the volume fraction of discrete 
particles and particle—particle collisions. Two-phase detonation is a highly transient 
two-phase flow problem involving strong shock waves; hence, thermal and mechan- 
ical non-equilibrium between the gas and particles should be considered. Further, 
uniform distribution of temperature is considered within particles due to their small 
particle sizes, and ideally, all reaction heat is absorbed by gas only. Based on the 
above assumptions, the two-dimensional governing equation of the gaseous phase in 
an Eulerian—Eulerian framework can be expressed as follows [9]: 


dU dF IG 
—+——+ —=S+W, 8.4 
T + F + T + (8.4) 


where U is the vector of conserved variables, F and G the conservation flux vectors 
in the x- and y-directions, S the vector of two-phase interaction source terms, and W 
the vector of chemical reaction source terms, respectively. They can be expressed as 


pı piu piv © 
U = Pns = Praf ‘ G =— PnsV : Ww = Wns ‘ 
pu pu’ +p puv 0 
pv puv pv’? +p 0 
E (E + p)u (E + p)v 0 
(8.5) 
and 
0 
Jp 
S= : , (8.6) 
0 
—fs+upJp 
—fy + VpIp 
=4p 7 (up fx + vp fy) + (Ep/Pp) Jp 
In the above equations, p; is the mass density of species i andi = 1, ..., ns; ns is 


the number of species contained in the gas mixture; p, u, and v are the gas pressure 
and x- and y-components of gas velocity, respectively; œ; is the mass production rate 
of gaseous species i by chemical reactions; and p and E are the mass density and 
total energy per unit volume of the gas mixture, expressed as 
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ns 


1 
= is E = h igi = 5 5 $ 8.7 
p 20e p P+ 5e(u +x?) (8.7) 
where A is the specific enthalpy of the gas mixture, calculated by 


SN h, (8.8) 


i=l 


with the specific enthalpy of each individual species h; as a function of gas tempera- 
ture T, obtained from a species thermodynamic data base, such as the NASA Glenn 
data base [10]. By further assuming each species performs as a perfect gas, the 
equation of state of the gas mixture is then given by 


ns 


Ro 
p= 2, ART, R= a, (8.9) 
where Ro = 8.314 J/(mol K) is the universal gas constant, and W; is the molar mass of 
species i. Variables in Eq. (8.6) are related to the properties of the discrete phase and 
the interaction between two phases. Their definitions will be given in the following 
paragraphs. 

As for the discrete phase, the governing equation has a similar form as Eq. (8.4): 


= Sp, (8.10) 


where U,, F, and G, are the vectors of conserved variables, conservation fluxes in 
the x- and y-directions of the discrete phase, which can be given by 


Pp PpUp PpYp 
PpUp Ppt, PpUpVp 
Up= | pvp |; Fp =] Ppüpyp |» Gp = Poy; , 
Ep Epup Epp 
Np NpUp NpVp 
-J, 
fc — UpIp 
Sp = fy —YpJp (8.11) 
Ip + (up fx T vp fy) = (Ep/Pp) Jp 
0 


In Eqs. (8.6) and (8.11), J, is the mass regression rate (the combustion rate) of the 
discrete phase and it is determined by the combustion model of the discrete phase. 
Pp is the density of the discrete phase in the suspension. up and v, are the x- and 
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y-components of the velocity of the discrete phase, respectively. E, and N, are the 
total energy per unit volume and the particle number density of the discrete phase, 
respectively. They can be calculated by 


6Pp 
TPmd? 


1 
= pp(us,+ v2), Np = 


: (8.12) 


Ep = Ppep + 


where e, is the specific internal energy of the discrete phase and can also be obtained 
from a species thermodynamic data base as a function of discrete phase temperature 
Tp; dp is the diameter of the particle; and p, is the material density of the discrete 
haan. 
Additionally, the x- and y-components of drag force acting on the discrete phase, 
fx and fy, can be modelled as follows: 


f= = Cod |V — V,|(u — up) 
; (8.13) 
i= Now Cpd;p|V — Vp|(v — vp) 
where Cp is the drag coefficient, 
24 1 p.2/3 
Cp = Re, (1 + sRe7/ ), for Re, < 1000 l (8.14) 
0.424, for Re, > 1000 


In Eqs. (8.13) and (8.14), the relative velocity and the relative Reynolds number 
Re, between the gas and discrete phases can be calculated by 


V-V, = w= ob v—-v aye (8.15) 
P r P 
and 
dp iag p 
Re, = E (8.16) 


where u is the viscosity coefficient of gas. Accordingly, the convection heat transfer 
between the two phases is expressed as follows: 


Gp = NprdpANu,(T — Tp), (8.17) 


with the two-phase Nusselt number expressed as functions of Re, and Prandtl number 
Pr: 


0.33 
Nu, = 2 + 0.459 Re? Pr . (8.18) 
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The CESE method has been proven to be of high accuracy and good stability in 
solving gaseous detonation problems in Sect. 8.1, and it is chosen to solve the two- 
phase detonation problems as well. That is, Eqs. (8.4) and (8.10) in the Eulerian- 
Eulerian framework can be solved using a CESE method, such as the a-a scheme. 
Theoretically, the source terms in the governing equations can be addressed together 
with the space-time integration in the CESE method [11-15]. However, a separated 
treatment of source terms is always employed in solving two-phase detonations, 
because the Jacobian matrixes of the interphase interaction and chemical reaction 
source terms in coupling treatment are rather complicated to calculate. Notably, 
with the source terms treated separately, it has been proven that the good accuracy 
and stability of the CESE method are preserved in high-speed reactive flow simula- 
tions [5, 16, 17]. On the other hand, the characteristic time scales of the interphase 
interaction and chemical reaction source terms are much smaller than that of flow 
dynamics, always leading to stiffness problems in two-phase detonation simulations. 
To overcome this problem, the operator-splitting technique with multiple sub-time 
steps [18] is always employed and then the source terms of interphase interactions 
and chemical reactions are explicitly integrated as ordinary differential equations. 
The detailed implementation process under the Eulerian—Eulerian framework can be 
illustrated as follows: 


CESE ~ — 
Un SS tat Ar = At/N 
y FBG Unt = Un 
=e +1 ( _T 
? Sp gi U, n+l — Up,n+1 


ua Uy > S®™ wm sm 

E 5 F F p 

=> 4 UT? = UM + Ar [S™ + WO] 
Ue = Ue + Ar S™ 


p.n+l p,n 
— FW) 
Unti = Una 8.19 
U —~y™ ’ (8.19) 
potl ~ “pnt 


where the subscripts n and m refer to the global time step (Az) and the sub-time step 
(Afr), respectively, and N is the total number of sub-time steps within one global 
convection time step of the CESE method. Depending on the degree of stiffness in 
the problem, N can be chosen to be 10—20. 

To test the accuracy and robustness of the CESE method in two-phase detonation 
simulations under the Eulerian—Eulerian framework, Wang et al. [15] applied the 
CESE method to simulate the detonation synthesis of titania (TiO2) nanoparticles. 
The corresponding chemical reaction can be described as follows: 


TiCl4(g) + O2(g) + 2H2(g) > TiO2(s) + 4HC1(g) (8.20) 
As seen, it involves TiO, solid particles as well as TiCl4, O2, Hz and HCI gases, 


implying that it is a gas—solid two-phase detonation problem. Simulation results 
showed that the simulated detonation profiles of gas density, pressure and temperature 
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Fig. 8.8 Profiles of gas density, velocity, pressure and temperature in the detonation synthesis of 
nanosized TiO? particles. Courtesy of G. Wang [15] 


are similar to those in the gaseous detonation wave, but a sudden drop of gas density 
can be captured after the detonation front because of phase transition from gas to 
solid particles, as shown in Fig. 8.8. Additionally, it is also found that the simulated 
particle size of the produced TiO, and the simulated peak pressure agreed well with 
other calculation and experiment data [19], respectively. 

As for gas-liquid two-phase detonation, Wang et al. [3] solved the two-phase 
planar detonations of liquid hydrocarbon fuels with the same CESE method under the 
Eulerian—Eulerian framework. It is also found that the detonation profiles are similar 
to those of gaseous detonations, but with higher gas density, pressure and temperature 
behind the detonation front due to higher energy density of liquid fuels. Additionally, 
the length of the reaction zone is found to be longer as well, which is relative to the 
slower combustion rate of the liquid particle. Table 8.3 shows the comparisons of 
simulated detonation speeds to the experimental values. It can be revealed that the 
relative errors are less than 10% for particle sizes small than about 145 um. As for 
two cases with larger particle sizes, the cases with relative errors exceeding 10%, it 
may be caused by the neglect of deformation of large liquid particles in modelling, 
since the combustion rate of liquid particle is influenced by particle deformation. All 
the above numerical results showed that that the CESE method is of high accuracy 
and good robustness in two-phase detonation simulations under an Eulerian—Eulerian 
framework. 

The Eulerian—Eulerian framework is a simple and effective way to deal with 
the discrete phase in two-phase detonations, but it has some inherent limitations in 
modelling realistic two-phase suspension in industries or experiments with a specific 
particle size distribution, where a relatively wide range of particle diameters are 
involved [20-22]. Under the Eulerian—Eulerian framework, all particles within one 
numerical mesh are assumed to be in the same states, such as the same particle size, 
temperature, velocity and so on. However, the number of particles within one mesh 
may be large, and their states may differ depending on their initial states and interac- 
tion histories with gas phase. Moreover, the forces and the heat transfers between the 
gas and particles also differ by the particle size, resulting in different temperatures and 


110 8 Application: Detonations 


Table 8.3 Comparisons of detonation speeds of liquid hydrocarbon fuels between simulations and 
experiments. Courtesy of G. Wang [3] 


Fuel Particle size Equivalence Experiment Simulation Relative error 
(um) ratio (m/s) (m/s) (%) 
CoHis | 20-30 0.41 1670 -7.50 
C6H14 20-30 0.49 1720 1671.30 —2.83 
C6H14 20-30 0.56 1700 1773.92 4.35 
CoHi4 20-30 0.68 1780 1935.89 8.76 
C10H20 145 1.0 2130 1946.93 —8.59 
Cj0H20 375 0.914 1810-1850 1996.72 7.93-10.32 
C10H20 1300 0.23 970-1250 1169.97 —6.40-20.62 


velocities of particles within one computational mesh. Further, it had been demon- 
strated that many two-phase detonation characteristics are significantly influenced 
by the particle size. Consequently, the Eulerian—Eulerian framework is insufficient to 
reflect the true physics of realistic two-phase detonations with particle size distribu- 
tions and to simulate them accurately, which raises the demand of modelling realistic 
two-phase detonation under the Eulerian—Lagrangian framework where the discrete 
particles are tracked individually. 

Under the Eulerian—Lagrangian framework, Eq. (8.4) is still applied as the 
governing equation of the gas phase, and U, F, G and W have the same forms 
as Eq. (8.5). However, the two-phase interaction source term S is expressed in the 
following form instead, 


0 
Tae 
1 D : 
S= — : , 
dV = 0 


— fk + U pk I pk 
— fyk + Vpk I pk 
— pk — (upk fsk + Vpk fyk) + tu? + E) -Jpk + epk ` Jpk 
(8.21) 


where subscript k represents all the quantities related to the kth particle (solid or 
liquid). To include all effects of particles into the source term S of the gaseous 
equation, the summation is done within the gaseous mesh element dV, and np is the 
number of particles in dV. 

The motion of every particle is then descripted using Newton’s laws of motion 
instead of the Eulerian form Eq. (8.10). For the kth particle, the corresponding 
governing equation can be written as 
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= Six, (8.22) 


where L,,; and Sp, are the vectors of the Lagrangian variables of the kth particle and 
the corresponding source terms, respectively, and they are expressed as 


T 
Lik = [Mm pk: Xpks Ypk, Mpk pk, MpkV pk, Ep] , (8.23) 
and 
T 
Spk = [- Jp, U pk, Vpk, Jrk, Sok —C pk J pk F apk] . (8.24) 
In Eq. (8.23), mpx and Ep; are the mass and total internal energy of the kth particle, 
respectively, and Epk = mpk-€pk- 
Further, the Lagrangian forms of the drag force and the convection heat of the 


kth particle can be easily derived from the according Eulerian forms Eqs. (8.13) and 
(8.17), as follows: 


fsk = = Cord2,p|V — Vpx|(u — upk) 


- , (8.25) 
Far = g Cond ph |V — Vor) (v — vor) 
and 
qpk = xdpkÀNup(T — Tpk). (8.26) 


Under the Eulerian—Lagrangian framework, the CESE method is again applied 
to solve the gas phase equation, while the source terms of interphase interactions 
and chemical reactions are treated separately and integrated explicitly as ordinary 
differential equations, along with the integration of particle Lagrangian equations, 
by using the operator-splitting technique as well, as depicted below, 


Ať = At/N 

CES as (0) T 

U, — U >4\U =U 
n S=W_o n+l Ne n+1 


Lakn = L pkn 


(m) (m) ( (m) 
Unti Dokar s™ wm, Spk 
=>) UST? = UM) + Ar[s® + wW] 
(m+1) _ 7 (m) (m) 
| rang = | ene + At’. Spk 
Un = UÑ 
sy tw (8.27) 
pkn+l = H pk,n+1 
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Notably, the Eulerian—Lagrangian framework had rarely been developed to solve 
high-speed reactive two-phase flows, mainly because the number of fine particles 
has proven too large and simulations too expensive to achieve in the past. Nowa- 
days, with the rapid development of computer technologies, parallel computation 
techniques such as the Message Passing Interface (MPI) technique may help to 
make simulations of two-phase detonations with fine particles under an Eulerian— 
Lagrangian framework possible. However, the application of MPI parallel compu- 
tation technique to an Eulerian—Lagrangian framework is not as straightforward as 
that under a purely Eulerian framework. A large amount of information exchange 
from one computational core to another is needed for calculations of the interaction 
source terms between two phases, leading to formidable communication cost even 
exceeding the computational cost, especially when the Eulerian framework is stag- 
gered with the Lagrangian framework due to the relative movement between these 
two phases. 

To solve the communication problem of parallel computing technique of high- 
speed two-phase flows under an Eulerian—Lagrangian framework, the traditional 
static data structure, such as multi-dimensional array, should be avoided using to store 
the information of the discrete phase described under the Lagrangian framework. 
Notably, the order of particles presented under the Lagrangian framework is not 
important and does not need to be preserved as its initial order. The only operation 
required for coding is to traverse every particle one by one. Inspired by the above 
facts, dynamic data structures can be introduced to store particle information and 
solve the communication problem under the Eulerian—Lagrangian framework. For 
example, the structures and the corresponding operations of linked lists used to store 
particle information are schematically shown in Fig. 8.9. As seen, each CPU owns 
one linked list to store information of the particles that are at the corresponding 
locations to the gas phase. In each linked list, one node represents one particle and 
consists of two parts: the data part that stores particle information, and the pointer 
that points to the next node of the particle for traversing. The pointer of the last node 
always points to the “NULL”, implying that the linked list has ended. 

Additionally, there are four basic operations for the linked list to adjust the 
particle’s storage location: allocate, free, delete and insert, as depicted in Fig. 8.9. 
The “allocate” operation is used to allocate new memory to store information about a 
“new” particle, the “free” operation to free the memory that stores information about 
an “old” particle, the “delete” operation to disconnect one particle from the linked 
list, and the “insert” operation to connect one particle at the end of the linked list. 
When the particle phase is staggered with the gas phase because of relative motion, 
the information of the particles, whose locations exceed the corresponding location 
ranges assigned to the present CPUs, will be transferred to and stored in the CPUs 
with the correct particle location ranges. For example, for CPU A in Fig. 8.9, where 
one particle is removed, the following sequence of operations is needed: 


(1) send the information of the specific particle to CPU B; 
(2) delete the particle node from the linked list; 
(3) free the memory of the separated node. 
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Fig. 8.9 Linked lists and operation sequences. Courtesy of Z. J. Zhang [9] 


Meanwhile, for CPU B in Fig. 8.9, where a new particle is received, the 
corresponding operation sequence is: 


(1) receive the information of the new particle from CPU A; 
(2) allocate a new node and store the information into the data field of the node; 
(3) insert the new particle node at the end of the linked list. 


With this dynamic data structure and the corresponding operation sequences, the 
information about Lagrangian particles is always stored in the CPUs of the correct 
Eulerian coordinates, and therefore excessive communication between CPUs when 
calculating gas-particle interactions is avoided, as shown in Fig. 8.10. Moreover, 
with the limitation of the global time step by the CFL condition, only “one” particle 
at most will cross the CPU boundary at every iteration; that is, the information of 
“one” particle at most will be transferred to the other CPU at one iteration step by the 
above operation sequence. As a result, the communication cost of the MPI parallel 
for the gas-particle interaction calculation will be reduced from O(N) to O(1) when 
using this data structure. Here, N is the number of particles stored in each CPU. 

To demonstrate the MPI parallelization performance with the use of the above 
linked lists and the corresponding operation sequences, Fig. 8.11 depicts the speedup 
parameters of a 2D Al-air detonation propagation problem using the CESE method 
under an Eulerian—Lagrangian framework, which involves approximately 24 million 
Eulerian meshes and about 75 million Lagrangian particles in the computational 
domain. This test case is simulated on the Tianhe-2 supercomputer from China with 
core numbers of 1, 2, 3, 4, 6, 12, 24, 48, 96, 192 and 384. The use of one core means 
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Fig. 8.10 Dynamic data structure using the linked list in the Eulerian—Lagrangian framework. 
Courtesy of Z. J. Zhang [9] 


that the simulation is done serially. As depicted in Fig. 8.11, when 384 cores are used, 
the code using linked lists still has a reasonable parallel efficiency of about 50% for the 
tested problem, while MPI parallelization under the Eulerian—Lagrangian framework 
using multi-dimensional static arrays is impossible, even with only 2 cores, as the 
communication cost is shown to be unacceptably large. This means that the linked list 
dynamic data structure works well in the MPI implement when solving two-phase 
detonations under the Eulerian—Lagrangian framework. 

To test the performances of CESE method in simulations of two-phase detonation 
under the Eulerian—Lagrangian framework, Shen et al. [23] simulated the gas—liquid 
two-phase detonations in the CioH22—07/air systems with different fuel droplet sizes 
and equivalence ratios. A deficit in the detonation speed compared to the corre- 
sponding purely gaseous one was observed in the gas-liquid suspensions with lean 
fuel and larger droplet sizes, while an increase in the detonation speed was observed 
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Fig. 8.12 Distribution of particle size near the detonation front in gas-liquid two-phase detonations. 
Courtesy of H. Shen [23] 


with very rich fuel. This special two-phase detonation feature in rich fuel is due to the 
slow combustion rate of the droplet, leading to a large amount of fuel burnt behind 
the C-J point, as shown in Fig. 8.12. As a result, the efficient equivalence ratio is 
reduced and detonation speed shifts to the equivalence ratio direction. The accuracy 
and robustness of the CESE method applied to two-phase detonation simulations 
under an Eulerian—Lagrangian framework had been demonstrated. 

Notably, simulations of the monodisperse two-phase detonation problems (with 
single particle size) under the Eulerian—Lagrangian framework always present the 
same results as those performed under the Eulerian—Eulerian framework. The raising 
of modelling two-phase detonations under the Eulerian—Lagrangian framework is in 
fact to simulate the realistic polydisperse two-phase detonations, where a particle 
size distribution is involved and the particle diameters are in a relatively wide range. 
One of the frequently used particle size distribution models to describe the polydis- 
perse suspension is the log-normal distribution function, expressed by the number 
frequency distribution function as follows [24]: 


1 1 /Ind, -Inday \? | 1 
n dp = P i 8.28 
f ( o) V2 00 z( 00 ) dp ( 


where d,,y and cog are the number median diameter and standard deviation of the 
distribution, respectively. In practice, the specific particle size distribution of a poly- 
disperse suspension is always given by other two particle size parameters that can be 
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easily measured, namely the volume-average diameter (d) and the mass-weighted- 
average particle diameter (dm). For example, the measured polydisperse parameters 
of the Al powder tested in the experiment of Zhang et al. [20] are d =2 um and dm 
= 3.3 um. The relationships between (dnm, o0) and (d, dm) can be obtained through 
the integrations of Eq. (8.28), as follows 


= œo 1/3 a 
d= p fas] = dye? 
0 


din = [f noa) f tacos] = me 
0 0 


Consequently, the corresponding parameters of the above-mentioned Al powder 
are dam = 1.37 wm and oy = 0.5. However, it is a good way to use the parameter set 
of (d, oo) to discuss the polydisperse suspensions in numerical simulations, because 
the effects of particle size distribution on polydisperse detonations can be easily 
discussed by varying o. Some specific particle size distributions with fixed d = 
2 um and different oo values are depicted in Fig. 8.13. 

In the followings, the typical Al-Air detonation problems with different particle 
size distributions, taking as an example, are simulated by the above numerical algo- 
rithm under the Eulerian—Lagrangian framework, to show the capacities of the CESE 
method in solving polydisperse two-phase detonation and show the differences 
between monodisperse and polydisperse two-phase detonations. The single-step 
global chemical reaction occurs in the Al-air suspension is 


(8.29) 


2Al(s) + 3/202(g) —> Al203(s) (8.30) 


f(d [X10-8/m] 


Fig. 8.13 Log-normal particle size distributions with d = 2 um. Courtesy of Z. J. Zhang [25] 
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To model the combustion rate of every single Al particle, the surface-kinetic- 
oxidation and diffusion hybrid combustion model (originally proposed by Zhang 
et al. [26]) can be employed. The Lagrangian form of the combustion rate of the kth 
Al particle is given as 


PT ae vaWar  kakksk 
pk,Al = pk, Al~O2 


: . 8.31 
"Vo, Wo, kak + ksk m 


where the reaction rates for the diffusion-controlled and kinetic-controlled combus- 
tion regime are expressed by 


vo, Wo,  Paidpk,Al 

2 
vaiWar 2Ciota K dko, al 
— Ea / RoTsk 


1/3 
n= (: + 0.276Re!/? Pr) 


í (8.32) 


ksk = koe 


In Eqs. (8.31) and (8.32), va; and vo, are the stoichiometric coefficients for Al and 
Op, respectively; Co, and Ciota are the mole concentrations of Oz and gas mixture, 
respectively; Ts is the particle surface temperature; and K, ko and Ea are model 
constants. 

As depicted in Fig. 8.14. Some special two-phase detonation features corre- 
sponding to multi-phase interaction can be clearly identified in the detonation front 
structures of monodisperse Al-air detonation. The first feature is known as “double 
peaks” in the gas pressure, density and velocity profiles (Fig. 8.14a, c), which is 
distinctly different from the single peak feature observed in gaseous detonations. 
It is demonstrated, from the one-dimensional flow theory in gas dynamics that the 
second peak in the detonation front structures of monodisperse Al suspension is 
caused by the dominant stage of heat transfer due to intense phase transition (Al 
evaporation) at a specific location after the shock front. The second feature, shown 
in Fig. 8.14b, is the plateau of particle temperature due to Al evaporation, which is 
equal to 2750 K and results in the observed “kink” in the gas temperature profile due 
to the intense heat transfer between the gas and particles. The third feature, shown 
in Fig. 8.14c, is the particle velocity lag in the velocity relaxation process, resulting 
in the alternative forces acting on particles and momentum transfers between the gas 
and particles. All these two-phase features for monodisperse Al-air detonation are 
similar with those obtained by Zhang et al. [26] and Teng and Jiang [27, 28]. 

However, as for the polydisperse detonation with a log-normal particle size distri- 
bution of oo = 0.5, most features of monodisperse Al detonation, including the 
double peaks in gas pressure, density, velocity profiles and the kink in gas tempera- 
ture profile, disappear, and only single peaks of gas quantities exist in the detonation 
front, which are quite similar to the wave front structures in gaseous detonations. The 
double peaks in the detonation front is demonstrated theoretically to be attributed to 
the space-dispersed phase transition processes of particles of different sizes result in 
an overall moderate heat transfer intensity, which hinders the formation of the heat- 
transfer-dominant stage. These differences between monodisperse and polydisperse 
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Fig. 8.14 Comparison of front structures in gas phase of polydisperse and monodisperse Al-air 
detonations: a pressure, b temperature, ¢ velocity, d density, e mass fraction of O2. Courtesy of Z. 
J. Zhang [25] 


two-phase detonations are relative to the multiple timescales and length scales in 
polydisperse suspension with a continuous particle size distribution. 

As for the multi-dimensional two-phase detonation features, Fig. 8.15 shows the 
comparisons of cellular detonation flow fields between monodisperse and polydis- 
perse detonations. As indicated in Fig. 8.15a, typical cellular detonation structures, 
including pairs of triple points, Mach stems, incident shock waves and pairs of trans- 
verse waves, can be observed in both the monodisperse and polydisperse detonation 
fronts. Moreover, the gas temperature and the Oz species mass fraction distributions 
for both monodisperse and polydisperse cases are non-uniform behind the detona- 
tion fronts, which are characterized by irregular local (high and low) temperature 
and O7 concentration zones, respectively, as shown in Fig. 8.15b, c. These irregular 
distributions of the flow field parameters in monodisperse and polydisperse detona- 
tions are both caused by the periodical motions of triple points along the detonation 
fronts. All these features in two-phase detonations are similar to those observed in 
gaseous detonations, except for the transverse waves. The transverse waves in both 
the monodisperse and polydisperse detonation fronts are weak and degenerate fairly 
fast in the rear flows, which are different from the strong transverse waves observed 
in purely gaseous detonations. According to Zhang et al. [26], these weak transverse 
waves can be attributed to the slow diffusion-controlled combustion of the majority 
of Al particles after their kinetic-inductions and a considerable amount of condensed 
Al oxide formed in the detonation products without contributing to gas pressure. 
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Fig. 8.15 Comparison of flow fields in Al-air detonation fronts. Left: monodisperse with dp = 
2 um at t = 0.6 ms. Right: polydisperse with d = 2 um and co = 0.5 at t = 1 ms. Courtesy of Z. 
J. Zhang [25] 


Further, larger detonation cell sizes are expected in polydisperse two-phase deto- 
nation compared to monodisperse detonation, since the reaction zone is larger, as 
shown in Fig. 8.16 by peak pressure contours. The estimated cell sizes of the monodis- 
perse detonation and the polydisperse detonations (oo = 0.5 and 0.8) are monodisperse 
= 10.5 + 0.5 mm, Aos = 13.3 + 0.8 mm and Apg = 20.0 + 1.8 mm, respectively. 
Ao.s is 27% larger than monodisperse, and Ao.g is even 190% larger. Figure 8.17 plots 
these detonation cell sizes as a function of the square of standard deviation og bya 
logarithmic scatter diagram, together with those of other three polydisperse detona- 
tions with different og. An approximately linear relationship between In A and og 
can be captured. Then, linear fitting is employed, and the fitting function appears to 
be Ind = 0.93670 + 0.0655, implying that the detonation cell size 4 of polydis- 
perse suspensions is an exponent function of the square of standard deviation of the 
distribution OG. This is an important quantitative relation of polydisperse two-phase 
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Fig. 8.16 Comparison of cellular detonations (peak pressure contours) of Al-air mixtures: a 
monodisperse with dp, = 2 um, b polydisperse with d = 2 ym and og = 0.5, and ¢ polydisperse 
with d = 2 um and oo = 0.8. Courtesy of Z. J. Zhang [25] 


Fig. 8.17 Detonation cell 2.5 
size à of Al-air suspension as 
a function of square of sis Gas 
standard deviation of, 2 with R? = 0.9929 


In 4 = 0.9367 of + 0.0655 


Ind 


1:5 


detonation to evaluate the effect of particle size distribution on detonation cell size. 
From the above results, the capacities of the CESE method in polydisperse two-phase 
detonation simulations under the Eulerian-Lagrangian framework is demonstrated 


as well. 
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Chapter 9 A) 
Other Applications get 


The CESE method has been applied to a wide range of scientific and engineering prob- 
lems since its inception in the 1990s. Although solving CFD problem is the primary 
goal of the CESE method, this general approach is actually applicable to a variety of 
PDE systems with physical backgrounds different from fluid dynamics. This chapter 
mainly introduces the application highlights of CESE in several representative fields. 


9.1 Hypersonic Aerodynamics 


Typical high-speed aerodynamic problems can be solved by using the CESE methods. 
Successful examples of such applications include the supersonic flow over a blunted 
flat plate [1], Mach reflection of shock waves [2], unsteady viscous flows in rocket 
nozzles [3], counterflow jets for the reduction of aerothermal load on spacecrafts [4, 
5], and the flow over roughness elements in a supersonic boundary layer [6]. 

With the increasing number of Mars exploration missions in recent years, the 
research interest of hypersonic aerothermodynamics has been revived. In a typical 
atmospheric-entry flight, the aircraft will encounter hypersonic gas flow, and the flow 
field is characterized by strong shock waves and thermochemical non-equilibrium 
effects. CFD has become a useful tool for in-depth understanding of complex 
flow physics and for complementary experimental prediction of aerodynamic 
characteristics. 

Using a 2D axisymmetric CESE solver based on hybrid meshes, Wen et al. [7— 
9] conducted a systematic study of hypersonic chemically reacting non-equilibrium 
flows over spheres. Three different working gases: nitrogen, air, and carbon dioxide 
were used in their simulations. The thermochemical non-equilibrium effects consid- 
ered therein include the vibrational energy relaxation of gas molecules as well as the 
dissociation and recombination reactions in the gas flows. The physically consistent 
coupled vibration—chemistry—vibration (CVCV) [10] two-temperature model was 
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employed. Here, the performance of the CESE method are shown by three different 
flow cases. The working gas, sphere radius, and the CESE simulation results [9] 
of the dimensionless shock stand-off distances (defined as the ratio of the stand-off 
distance to the sphere radius) of each case are listed in Table 9.1, along with the 
experimental and theoretical results of Wen and Hornung [11]. The CESE results are 
shown in good agreement with experimental and theoretical results. Further investi- 
gation of the shape of shock waves is presented in Fig. 9.1, in which the numerical 
density contours are overlaid on the experimental finite fringe differential interfer- 
ograms. The shapes of the simulated bow shocks with the CESE method match the 
experimental results well. 

Another important thermochemical non-equilibrium process in an atmospheric- 
entry flow is the ionization in a high-temperature shock layer. Massimi et al. [12] 
conducted a numerical investigation of hypersonic ionized air flows over rounded 
nose geometries, using a 2D axisymmetric CESE solver based on hybrid meshes. The 
multi-species NS equations were solved with the two-temperature seven-species Park 
model for ionization reactions. The weakly ionized air flow in the Radio Attenuation 
Measurement (RAM-C II) flight test at 71 km altitude was simulated and compared 
with experiments and numerical results of Candler and MacCormack [13]. The CESE 
results [12] of the maximum electron number densities along the direction normal to 
the vehicle’s surface agree well with the flight test measurements and the numerical 
results of Candler and MacCormack [13]. 


Table 9.1 The working gases, sphere radii, and dimensionless shock stand-off distances (last three 
columns) in different flow cases 


Case no Gas Radius (inch) Experiment [11] Theory [11] Simulation [9] 
1 N2 1.5 0.100 0.095 0.100 
2 Air 1.0 0.105 0.093 0.095 
3 CO2 2.0 0.088 0.084 0.087 


(a) Case 1 (b) Case 2 (c) Case 3 


Fig. 9.1 Comparison of experimental shapes of the bow shocks with the numerical predictions by 
CESE. Courtesy of C. Y. Wen [9] 


9.2 Aeroacoustics 125 


9.2 Aeroacoustics 


Thanks to the rapid development of CFD, it is now expected to solve aeroacoustics 
problems through high-resolution numerical simulations of unsteady compressible 
Euler or NS equations. However, in many computational aeroacoustic problems, 
the coexistence of shock waves and small disturbances (acoustic waves) imposes 
stringent requirements on numerical methods. 

The sound-shock interaction problem was used by Wang et al. [14] to examine 
the accuracy of the CESE method for aeroacoustic problems involving shock waves. 
Loh et al. [15] further considered three selected problems, namely a linear bench- 
mark problem, instability waves on a free shear layer, and shock—vortex interactions, 
to illustrate the feasibility of using CESE as a tool for computational aeroacoustics 
(CAA). Good numerical results with high resolution and low dispersion were demon- 
strated by their second-order CESE scheme. Additionally, the authors also pointed 
out that the non-reflecting boundary condition, which plays an important role in 
CAA, is much simpler to implement in the CESE method than in traditional methods 
(see [15, 16]). To deal with practical multi-dimensional CAA problems where large 
disparity in grid cell size is inevitable, Yen et al. [17, 18] modified the Courant number 
insensitive (CNI)-CESE scheme [17] and the local time-stepping CESE scheme [18] 
to improve the numerical efficiency. A 3D simulation was performed for the propaga- 
tion of a Gaussian acoustic pulse and a vorticity wave embedded in a Mach 0.5 mean 
flow. Once again, the CESE results were in good agreement with the corresponding 
analytical solution in preserving both the form and amplitude of the waves [17]. Yen 
[18] further conducted a series of 3D CESE simulations for the Helmholtz resonator 
problem. By using a second-order CESE scheme, excellent agreement between the 
linear acoustic theory and the CESE solution was achieved. 

The CESE method has also been employed as a numerical tool for CAA prob- 
lems in practical engineering applications related to engines and aircrafts. To shed 
light upon an aeroacoustic resonance phenomenon often encountered by conver- 
gent—divergent nozzles under transonic conditions, Loh and Zaman [19, 20] devel- 
oped an axisymmetric NS solver using the CESE method and conducted a numerical 
investigation. A 3D CESE NS solver was also implemented by Loh et al. [15] and 
applied to compute the screech noise generated by an under-expanded supersonic 
jet. Kim et al. [21] simulated the supersonic unsteady flow over an open cavity to 
investigate the mixing-enhancement and flame-holding capability in the scramjet 
engine. The CESE method solving 2D NS equations successfully captured the self- 
sustained oscillations in the supersonic cavity flows. The computed frequencies and 
amplitudes of the pressure oscillations compare favourably with the theoretical and 
experimental data. Cheng et al. [22] also studied subsonic and supersonic flows over 
open-cavity geometries by an unsteady NS solver based on the CESE method. When 
compared with experimental and analytical data, the generation and propagation of 
flow-induced acoustic waves were faithfully captured, and satisfactory results were 
obtained for the oscillation frequency of the dominant mode. 
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9.3 Solid Dynamics 


In recent years, high-resolution discontinuity-capturing numerical methods origi- 
nally devised for CFD problems have been utilized to solve nonlinear solid-dynamic 
problems with Eulerian formulations. Among these methods, the CESE method has 
attracted considerable attention from researchers in computational solid dynamics 
because of its simple logic, high accuracy, and ability to capture the behaviours of 
nonlinear waves. 

Wang et al. [23] investigated numerically two high-velocity impact problems in 
elastic—plastic materials: (1) Taylor copper-bar-impact problem and (2) the pene- 
tration of a long rod (tungsten heavy alloy) into a steel target. In this study, the 
CESE method combined with the level-set technique was adopted to solve Eulerian 
governing equations that describe the solid dynamics and to trace material interfaces. 
Chen et al. [24] extended the CESE solver developed by Wang et al. [23] to simu- 
late high-velocity impact problems involving elastic-plastic flows, high strain rates, 
and spall fractures. Excellent agreement was observed between their simulation of 
aluminium plates colliding with stainless-steel plates and the experimental data. 

As known, the governing equations of elastic waves in solids can be properly 
written as a set of fully coupled first-order hyperbolic PDEs, including the conser- 
vation laws of mass and momentum as well as the rate-type constitutive relations 
for materials, by treating the density, velocity, and stress components as primitive 
unknowns. Therefore, the CESE method, which is suitable for hyperbolic systems, 
can be intuitively applied to simulate linear and nonlinear waves in elastic solids. The 
CESE simulations of resonant standing waves arising from a time-harmonic external 
axial load and compression waves arising from a bi-material collinear impact was 
demonstrated by Yu et al. [25]. Another study on the propagation and reflection of 
extensional waves in an abruptly stopped elastic rod using the CESE method was 
carried out by Yang et al. [26]. Moreover, Chen et al. [27] performed CESE simula- 
tions to study planar-wave expansion from a point source in an anisotropic solid of 
cubic symmetry (gallium arsenide) and Yang et al. [28] extended the CESE simula- 
tions to stress waves in solids of hexagonal symmetry, illustrating wave propagation 
in a heterogeneous solid composed of three blocks of beryl with different lattice 
orientations. Lately, Lowe et al. [29] conducted a comprehensive study of nonlinear 
longitudinal waves (including both weak and strong shocks, rarefactions, and contact 
discontinuities) in tapered elastic rods using the CESE method as a numerical tool. 
In the above studies, the CESE simulations effectively captured waves in a variety of 
solid materials and provided results consistent with the available analytical solutions. 
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9.4 Magnetohydrodynamics 


The plasma flows in aerospace applications and astrophysics can be described and 
solved generally by the magnetohydrodynamic (MHD) equations, which combine 
the NS equations for fluid dynamics and the Maxwell equations for electromagnetics. 
Due to the divergence-free constraint usually imposed for the magnetic field B in 
the computational MHD problems, some special treatments are added into the MHD 
numerical schemes to enforce V - B = 0. Since the proposal of the CESE method, 
researchers have been attempting to solve MHD problems accurately and efficiently 
with the CESE method. 

Zhang et al. [30, 31] used the CESE method to study MHD benchmark prob- 
lems, including an MHD shock-tube problem and an MHD vortex problem. Their 
CESE results without additional treatments for V - B = 0 compared favourably with 
previously reported reference solutions. In fact, Zhang et al. [30, 31] performed the 
simulations with the baseline CESE scheme and the CESE scheme in conjunction 
with a special treatment to maintain V - B = 0. Nevertheless, no obvious difference 
was observed. 

Feng et al. [32] developed a numerical platform based on the CESE method to 
investigate solar—interplanetary physics and space weather. Three dimensional MHD 
equations were solved with the adaptive mesh refinement (AMR) to better resolve 
flow features that have spatial scales many orders of magnitude smaller than the vast 
size of solar—interplanetary space. This CESE-based MHD-simulation approach was 
validated through the numerical study of solar corona and solar wind and comparison 
with observation data. Recently, Yang et al. [33, 34] extended this CESE-based MHD 
solver to a high-order version [33] and an upwind version [34] based on the CESE 
modifications of Shen et al. [35] and Shen and Wen [36], respectively. Numerical 
tests of MHD-vortex and MHD-blast-wave problems demonstrated the accuracy 
improvement of the CESE results by applying both the high-order and the upwind 
extension. Furthermore, a new strategy to keep the magnetic field fundamentally 
divergence-free was proposed by Yang et al. [33], taking the advantage of the CESE 
method that the spatial derivatives of the magnetic field are treated as marching 
variables in the algorithm. 

Notably, the successful combination of MHD modelling and CESE simulation 
was elaborated in a recent monograph of Feng [37]. For more details about the 
CESE method for MHD, readers are recommended to refer to Feng [37] and the 
references therein. 
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Chapter 10 A) 
Summary get 


In this book, the space-time CESE method is introduced as a novel and distinctive 
numerical approach to solving conservation equations in engineering sciences. Also, 
from a historical perspective, the development of the CESE method since the 1990s 
has been reviewed and the remarkable improvements and extensions of the CESE 
method in the past few years has been summarized. Various applications of the CESE 
method have been presented to emphasize its numerical performance under different 
scenarios, including many research topics in engineering, such as compressible multi- 
fluid flows and detonations. 

The most elegant part of the CESE method might be the non-dissipative a scheme 
described in Chap. 2, where the unknown variable u and its spatial derivative uy 
can be obtained directly from the conservation laws. However, for the purpose of 
shock capturing, two families of CESE schemes have been developed on the basis 
of a scheme. The first family is called the central CESE scheme, including the a—a 
and CNI schemes. In these schemes, the spatial derivatives are updated by specially 
designed procedures, where artificial dissipation can be added. The central CESE 
schemes avoid any characteristic-based techniques (i.e., they are free of Riemann 
solver). The second family is the upwind CESE schemes, in which the upwind flux 
solver (e.g., Riemann solver) can be utilized, but it just plays a supporting role. 
By combining the ideas in the original CESE method and the conventional upwind 
FVM, the upwind CESE scheme managed to retain the forms of discretized equations 
for u and ux in the a scheme, while the numerical dissipation can be introduced in 
a reasonable manner. It is worth noting that the dissipation of the upwind CESE 
scheme is naturally insensitive to the CFL number. 

Two important issues about the CESE method, namely the CESE schemes on 
unstructured meshes and the high-order CESE schemes, have been discussed in 
Chaps. 4 and 5, respectively. The CESE method is naturally multi-dimensional and 
compatible with unstructured meshes. A high-order extension of the CE/SE method 
in both space and time can be achieved by storing and updating high-order derivatives 
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(€.g., Uxx and u,x,) at each solution point, without enlarging the stencil or adding stages 
in time integration. 

The CESE method has low dissipation and high compactness. When applied to 
the simulations of complex physical processes, the CESE method can catch shock 
waves, contact discontinuities, fine structures, and small disturbances, with high reso- 
lution and strong robustness. Therefore, the CESE method demonstrates good perfor- 
mances in the numerical simulations of wave-propagation problems (e.g., shock 
waves, acoustic waves, detonation waves, stress waves in solid, and the electromag- 
netic waves), interfacial instabilities, and shock—bubble/droplet interactions. In many 
hot research areas including hypersonic aerodynamics, shock dynamics, detonation, 
aeroacoustics, MHD, and solid dynamics, the CESE method proves to be suitable 
and shows a good development prospect. 

In recent years, an important topic in CESE research has been its application 
to highly nonuniform and high-aspect-ratio computational meshes, which are often 
required in simulations of boundary layers. In particular, the CESE method has been 
recently applied to aerodynamic heating problems and direct numerical simulation of 
laminar or turbulent flows. It is also worth noting that the CESE schemes are mainly 
applied to unsteady problems due to their explicit nature in time with the associated 
restriction on time step size. In order to solve steady problems efficiently, the exten- 
sions to implicit time-stepping CESE schemes are under developing. Further analyses 
of the numerical properties of CESE schemes will be very useful to solidify their 
mathematical foundation. The modified equation analysis and the spectral property 
analysis for the upwind CESE scheme are in progress now. 

Along with the continuous improvement of high-performance computing, numer- 
ical methods emerge one after another. Many challenging numerical simulations can 
be accomplished today. However, none of the existing numerical methods can be 
accurate, robust, and efficient under all situations. In this sense, the CESE method may 
not be superior to the existing methods, but it does provide an alternative approach 
which deserves more evaluations and further investigations. The authors wish this 
book can familiarize the readers with the CESE method. 
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Appendix 


1. The a-a CESE solver for solving the 1D scalar convection problem. 


double avg(double x1,double x2) 


{ 
double x, aa; 


aa=2.0;//6*fabs (x1-x2) / (fabs (x1)+fabs(x2)+1le-10); 


x= (pow (fabs (x1) ,aa)*x2+pow(fabs(x2),aa)*x1)/ 


(pow (fabs (x1) ,aa)+pow(fabs (x2) ,aa)+1le-20); 


return x; 
}// Compute weighted average of x1 and x2. 


double CFL (double U[N+3] ) 
{ 
int i; 
double maxvel,vel; 
maxvel=1e-10; 
for (i=1;i<=N+1;i++) 
{ 
vel=1.0; 
if (vel>maxvel)maxvel=vel; 
} 
return Sf*dx/maxvel; 
}// Compute timestep. 
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double const dx= (double (L2-L1) ) /double(N) ; //Mesh size 


include <stdio.h> 
include <math.h> 
include <stdlib.h> 
include <time.h> 
define end_time 2.0 // Computing time. 
define L1-1.0 // Left end of the domain. 
define L2 1.0 // Right end of the domain. 
define N 200 // Mesh number. 
define Sf 0.8 // CFL number. 
define MAX(x,y) (((x)>(y))?(x):(y)) 
define MIN(x,y) (((x)<(y))?(x):(y)) 
double U1 [N+3] , U2 [N+3] ,Ux1[N+3] ,Ux2 [N+3],Ut[N+3],F[N+3],Ft[N+3]; 
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void Bound (double U[N+3] , double Ux[N+3]) 
{ 

U[0]=U[N]; 

U[N+2]=U[2]; 

Ux[0]=Ux[N]; 

Ux[N+2]=Ux[2]; 
}// Apply boundary condition. 


void U2Fu (double U[N+3], double Fu[N+3] ) 


{ 
intea4 
for (1=0;1<N+3;1i++) 
{ 
Fu[i]=U[i]; 
} 


}// Compute flux. 


void ComputeTemporalDerivatives (double U[N+3] , double Ux[N+3]) 
{ 
Int i; 
for (i=0;i<=N+2; i++) 
{ 
Ut[i]=-Ux[i]; 
Ft[i]=Ut[i]; 
} 
}// Compute the temporal derivatives of U and F. 


void CESE_TimeMarch (double U[N+3] , double Ux[N+3] , double 
Fu[N+3] , double 
Fut [N+3] , double Unew[N+3] , double Uxnew[N+3] , double dt, int ishalf) 
{ 

IRELE? 

double Uxplus,Uxminus,Uleft,Uright,Fleft,Fright; 

for (i=1-ishalf;i<=N+1;i++) 

{ 

I=i+ishalf; 


Uleft=U[I-1]+dx*Ux[I-1]/4; 
Uright=U[I]-dx*Ux[I]/4; 
Fleft=Fu[I-1]+dt*Fut[I-1]/4; 
Fright=Fu[I]+dt*Fut[I]/4; 


Unew[i]=0.5*(Uleft+Uright+dt/dx*(Fleft-Fright)); 
Uxminus=2* (Unew[i]-U[I-1]-dt*Ut [I-1] /2) /dx; 
Uxplus=2* (U[I]+dt*Ut [I] /2-Unew[i] ) /dx; 
Uxnew[i]=avg(Uxminus,Uxplus) ; 


} 
}// Scheme marching from integer-step to half-step, or from half- 
step to integer-step. 


void Initialize (double U[N+3] ,double Ux[N+3] ) 
{ 

int i; 

for (i=0;i<=N+2; i++) 


{ 
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if (i > =(N/4+1)&&i<=(3*N/4+1)) 


{ 

U[i]=1; 
} 
else 
{ 

U[i]=0; 
} 
Ux[i]=0; 


} 
}// Initialization. 


void Output (double U[N+3],double t) 
{ 
int i; 
double x; 
FILE *fp; 
fp=fopen("result.txt", "w+"); 
for (i=1;i<=N+1;1++) 


{ 

x=L1+ (i-1)*dx; 

fprintf (fp, "%20.5e\t%20.5e\n",x,U[i]); 
} 
fclose(fp); 


}// Output result. 


void CESE_Solver () 


{ 
int Nstep=0; 
double dt,t=0; 
Initialize(U1,Ux1); 
while (t<end_time &&Nstep>=0) 
{ 
dt=CFL (U1); 
t=t+dt; 
Nstep++; 
printf ("Nstep=%d,dt=%e,t=%e\n",Nstep,dt,t); 
U2Fu (U1,F); 
ComputeTemporalDerivatives (U1,Ux1); 
CESE_TimeMarch (U1,Ux1,F,Ft,U2,Ux2,dt,1);//the 1st half 
time step 
U2Fu(U2,F); 
ComputeTemporalDerivatives (U2,Ux2) ; 
CESE_TimeMarch (U2,Ux2,F,Ft,U1,Ux1,dt,0);//the 2nd half 
time step 
Bound (U1, Ux1) ; 
} 


Output (U1,t); 
}// a-a CESE solver. 


int main(void) 


{ 


CESE_Solver(); 
return 1; 


136 


Appendix 


2. The a-a CESE solver for solving the Sod-shock problem. 


include <stdio.h> 
#include <math.h> 
include <stdlib.h> 
include <time.h> 
define end_time 0.1 
define L1 0.0 
define L2 2.0 
define LO 1.0 
define N 200 

define Sf 0.8 
define GAMA 1.4 
define roul 1.0 
define ul 0.0 
define p1 1.0 
define rou2 0.125 
define u2 0.0 
define p2 0.1 
define MAX(x,y) (((x 
define MIN(x,y) (((x 
double U1 [N+1] [3] ,U2 
Ux2[N+1] [3], Ft[N+1] 


double CFL (double U[N 


double const dx= (double (L2-L1))/double (N); 


// Computing time. 
// Left end of the domain. 
// Right end of the domain. 


// Mesh number. 
// CFL number. 
// Ratio of specific heats. 


// Initial conditions 
) > (y))? (x): (y)) 
)<(y)) ? (x): (y)) 
N+1] [3], Ut[N+1] [3],F[N+1] [3] ,Ux1[N+1] [3], 
lh 


) 
R 


// Mesh size. 


+1) [3]) 


{ 
int i; 
double maxvel=le-4,vel,p,u; 
for (i=1;1<=N-1;1i++) 


{ 
u=U[i] [1]/U[1i] [0]; 
p= (GAMA-1) * (U[i] [2]-0.5*U[i] [0] *u*u); 
vel=sqrt (GAMA*p/U[i] [0])+fabs(u) ; 
if (vel > maxvel)maxvel=vel; 
} 


return Sf*dx/maxvel; 
}// Compute timestep. 


double avg (double x1, double x2) 
{ 
double x,aa=1.0; 
x= (pow (fabs (x1) ,aa) *x2+pow (fabs (x2) 
(pow (fabs (x1) ,aa)=+pow (fabs (x2), 
return x; 
}// Compute weighted average of x1 and x2. 


,aa) *x1)/ 
aa)+le-20); 


void Init () 
{ 
int i,k; 
double x; 
for (1=0;i<=N;i++) 
{ 


x=L1+i*dx; 
if (x<=L0) 
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U1 [I] [0]=roul; 
U1[I][1]=roul*u1; 
U1[I] =p1/ (GAMA-1)+roul*ul*ul/2; 


N 


U1[I][0]=rou2; 
U1[I] =rou2*u2; 
UL[T] =p2/ (GAMA-1) +rou2*u2*u2/2; 


H 


N 


} 
for (i=0;i<=N; i++) 
{ 
for (k=0;k<3;k++) 
{ 
Ux1[i] [k]=0; 
} 
} 
}// Initialization. 


void Bound (double U[N+1] [3]) 
{ 
for(int k=0;k<3;k++) 


else 

{ 
U[N] [kK] =U[N-1] [k]; 
Ux1[N] [k]=0; 
U[0] [k]=U[1] [k]; 
Ux1[0] [k]=0; 


} 
} 
}// Apply boundary condition. 


void U2F (double U[N+1] [3], double F[N+1] [3]) 
{ 

int i 

double u,p; 

for (1=0;i<N+1;i++) 


{ 


u=U[i] [1] /U[i] [0]; 
p=(GAMA-1)*(U[i][2]-0.5*U[i][1]*U[i][1]/U[i][0]); 
F[i][0]=U[i][1]; 

F[i][1]=U[i] [0] *u*utp; 

Fi] [2]=(U[i] [2] +p) *u; 


} 
}// Compute flux. 


138 Appendix 


void ComputeDerivative (double U[N+1] [3] ,double Ux[N+1] [3] ) 
{ 

int i; 

double u; 

for (i=0;1<=N;1i++) 


=U[i] [1] /u[i] [0]; 

Ut [i] [0]=-Ux[i] [1]; 

Ut [i] [1] =- (GAMA-3) *u*u*Ux[i] [0] /2- (3-GAMA) *u*Ux[i] [1] 
- (GAMA-1) *Ux[i] [2]; 

Ut [i] [2] =- ( (GAMA-1) *u*u*u-GAMA*u*U[i] [2] /U[i] [0])* 
Ux[i] [0] - (GAMA*U[i] [2] /U[i] [0]-3* (GAMA-1) /2*u*u) *Ux[i] [1] - 
GAMA*u*Ux[i] [2]; 

Ft[i] [0]=Ut[i] [1]; 

Ft [i] [1] =(GAMA-3) *u*u*Ut [i] [0] /2+(3-GAMA) *u*Ut [i] [1] 
+ (GAMA-1) *Ut [i] [2]; 

Ft[i] [2]=( (GAMA-1) *u*u*u-GAMA*u*U[i] [2] /U[i] [0])* 
Ut [i] [0]+(GAMA*U[i] [2] /U[i] [0]-3* (GAMA-1) /2*u*u) *Ut[i] [1] 


+GAMA*u*Ut [i] [2]; 
} 
}// Compute the derivatives of F. 


void CESE_TimeMarching(double U[N+1][3], double Ux[N+1][3], 
double Fu[N+1] [3],double Fut [N+1] [3], double Unew[N+1] [3], double 
Uxnew[N+1] [3], double dt, int ishalf) 

{ 


int.i).j4 5; 
double Uleft,Uright,Fleft,Fright,Ux_minus,Ux_plus; 
for (i=1-ishalf;i<N;i++) 
{ 
T=itishalf; 
for (j=0;5<3;j++) 
{ 
Uleft=U[I-1] [j]+dx*Ux[I-1] [51/4; 
Uright=U[I] [j]-dx*Ux[I][j3]/4; 
Fleft=F[I-1] [j]+dt*Ft[I-1][j]/4; 
Fright=F [I] [j]+dt*Ft[1I][j1/4; 
Unew[i] [j]=0.5* (Uleft+Uright+dt/dx* (Fleft 
-Fright)); 
Ux_minus=2* (Unew[i] [j]-U[I-1] [j]-dt* 
Ut[I-1] [3]/2) /dx; 
Ux_plus=2* (U[I] [j]+dt*Ut[T] [j]/2-Unew[i] [j]) /dx; 
Uxnew[i] [j]=avg(Ux_minus,Ux_plus) ; 
} 
} 
}// Scheme marching from integer-step to half-step, or from half- 
step to integer-step. 


void CESE_1D_Solver () 
{ 


int i1,comput_num; 
double dt,t=0; 
FILE *fp; 
comput_num=0 ; 
Init(); 
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while (t<end_time) 
{ 
comput_num++; 
dt=CFL (U1); 
t=t+dt; 
UZE(UL, £)% 
ComputeDerivative (U1,Ux1); 
CESE_TimeMarching (U1,Ux1,F,Ft,U2,Ux2,dt,1); 
UZE (U2); 
ComputeDerivative (U2,Ux2); 
CESE_TimeMarching (U2,Ux2,F,Ft,U1,Ux1,dt,0); 


Bound (U1) ; 


} 
}// a-a CESE solver. 


void Output (double U[N+1] [3]) 
{ 
EN, 39 
double x; 
FILE *fp; 
double rou,u,p; 
fp=fopen("result.txt", "wt"); 
for (i=0;i<N+1;1i++) 


{ 


x=L1+i*dx; 
rou=U[i] [0]; 
u=U[i] [1]/U[i] [0]; 
p= (GAMA-1) * (U[i] [2]-0.5*U[i] [0] *u*u); 
fprintf (fp, "%20.5e\t%20.5e\t%20.5e\t%20.5e\n", 
x,rou,u,p); 
} 
fclose(fp); 
}// Output result. 


int main () 

{ 
CESE_1D_Solver(); 
Output (U1) ; 
return 1; 


